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

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

MATLAB 2016a Release Summary for Scientific Computing

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

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

Julia iFEM 2: Optimizing Stiffness Matrix Assembly

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 ith 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 ... READ MORE

Optimizing .*: Details of Vectorization and Metaprogramming

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


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 ... READ MORE

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);
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]);
using PyPlot;
fig = figure(figsize=(20,10));

This plots the logistic map. If you take the same code and change the array dereferencing from A[i,j] to A(i,j) you get MATLAB code. However, I couldn't get MATLAB to have it looks as good, but that's beside the point.

Julia iFEM1: Porting Mesh Generation

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