A Comparison Between Differential Equation Solver Suites In MATLAB, R, Julia, Python, C, Mathematica, Maple, and Fortran

September 26 2017 in C, Differential Equations, Julia, Mathematica, Mathematics, MATLAB, Programming, R | Tags: | Author: Christopher Rackauckas

Many times a scientist is choosing a programming language or a software for a specific purpose. For the field of scientific computing, the methods for solving differential equations are one of the important areas. What I would like to do is take the time to compare and contrast between the most popular offerings. This is a good way to reflect upon what’s available and find out where there is room for improvement. I hope that by giving you the details for how each suite was put together (and the “why”, as gathered from software publications) you can come to your own conclusion as to which suites are right for you.

(Full disclosure, I am the lead developer of DifferentialEquations.jl. You will see at the end that DifferentialEquations.jl does offer pretty much everything from the other suite combined, but that’s no accident: … READ MORE

50 Comments

MATLAB 2016a Release Summary for Scientific Computing

March 15 2016 in MATLAB, Programming | Tags: , , , | Author: Christopher Rackauckas

There is a lot to read every time MATLAB releases a new version. Here is a summary of what has changed in 2016a from the eyes of someone doing HPC/Scientific Computing/Numerical Analysis. This means I will leave off a lot, and you should check it out yourself but if you’re using MATLAB for science then this may cover most of the things you care about.

  1. Support for sparse matrices on the GPU. A nice addition is sprand and pcg (Preconditioned Conjugate Gradient solvers) for sprase GPU matrices.
  2. One other big change in the parallel computing toolbox is you can now set nonlinear solvers to estimate gradients and Jacobians in parallel. This should be a nice boost to the MATLAB optimization toolbox.
  3. In the statistics and machine learning toolbox, they added some algorithms for high dimensional data and now let you run kmeans … READ MORE
Leave a Comment

Julia iFEM3: Solving the Poisson Equation via MATLAB Interfacing

January 24 2016 in FEM, Julia, MATLAB, Programming | Tags: , , | Author: Christopher Rackauckas

This is the third part in the series for building a finite element method solver in Julia. Last time we used our mesh generation tools to assemble the stiffness matrix. The details for what will be outlined here can be found in this document. I do not want to dwell too much on the actual code details since they are quite nicely spelled out there, so instead I will focus on the porting of the code. The full code is at the bottom for reference.

The Ups, Downs, and Remedies to Math in Julia

At this point I have been coding in Julia for over a week and have been loving it. I come into each new function knowing that if I just change the array dereferencing from () to [] and wrap vec() calls around vectors being used as … READ MORE

Leave a Comment

Julia iFEM 2: Optimizing Stiffness Matrix Assembly

January 23 2016 in FEM, Julia, MATLAB, Programming | Tags: , , , | Author: Christopher Rackauckas

This is the second post looking at building a finite element method solver in Julia. The first post was about mesh generation and language bindings. In this post we are going to focus on performance. We start with the command from the previous post:

node,elem = squaremesh([0 1 0 1],.01)

which generates an array elem where each row holds the reference indices to the 3 points which form a triangle (element). The actual locations of these points are in the array node, and so node(1) gives the points in the (x,y)-plane for the $$i$$th point. What the call is saying is that these are generated for the unit square with mesh-size .01, meaning we have 10201 triangles.

The approach to building the stiffness matrix for the … READ MORE

Leave a Comment

Optimizing .*: Details of Vectorization and Metaprogramming

January 21 2016 in Julia, MATLAB | Tags: , , , , , | Author: Christopher Rackauckas

For many of us mathematicians, we were taught to use MATLAB, and we were taught to vectorize everything. I mean obviously if we have matrices $$A$$, $$B$$, and $$C$$ and want to multiply element-wise (say to solve a reaction-equation at each point in space), then the optimal code is

A.*B.*C

No questions to ask, right? Actually, this code isn’t as optimized as you’d think. Lets dig deeper.

BLAS, Linpack, and MKL

The reason you are always told by “the lords of numerical math” to vectorize your code is because very smart programmers worked really hard on making basic things work well. Most of the “standardized” vectorized computations are calling subroutines from packages known as BLAS and LINPACK. To see what version your MATLAB is using, you can call

READ MORE

6 Comments

Quick Optimizations in Julia for Performance: A Practical Example

January 19 2016 in Julia, MATLAB | Tags: , , , , | Author: Christopher Rackauckas

Let’s take a program which plots the standard logistic map:

r = 2.9:.00005:4;
numAttract = 100;
steady = ones(length(r),1)*.25;
for i=1:300 ## Get to steady state
  steady = r.*steady.*(1-steady);
end
x = zeros(length(steady),numAttract);
x[:,1] = steady;
for i=2:numAttract ## Now grab some values at the attractor
  x[:,i] = r.*x[:,i-1].*(1-x[:,i-1]);
end
using PyPlot;
fig = figure(figsize=(20,10));
plot(collect(r),x,"b.",markersize=.06)
savefig("plot.png",dpi=300);

This plots the logistic map. If you take the same code and change the array … READ MORE

8 Comments

Julia iFEM1: Porting Mesh Generation

January 19 2016 in FEM, Julia, MATLAB, Programming | Tags: , , , , | Author: Christopher Rackauckas

My first project on the quest for a Julia finite element method is a simple homework problem. Just some background, this is for UC Irvine’s graduate Computational PDEs 226B course where in the first quarter we did all sorts of finite difference methods and now is our first foray into finite element methods. The purpose of the project is to grasp the data structure enough to use simple tools (i.e. mesh creation and plotting) to create a finite element solver for Poisson’s equation in 2D and check the performance differences. For those who haven’t programmed much this is a great learning exercise, but being a pretty standard exercise the MATLAB code took no time to writeup and so I started thinking: how does Julia compare to MATLAB for solving simple PDEs with the finite element method?

Testing this would take … READ MORE

Leave a Comment