The choice - 1 (Intro)

July 22, 2011 at 03:07 - by eelvex   |   one comment (Leave a comment)

Share on RedditShare on FacebookShare on Google+Tweet about this on TwitterPrint this page

[ Get the code at github ]

A programming language is low level when its programs require attention to the irrelevant.

Alan J. Perlis

It's not only that physicists can't program (we rarely use version control, we don't plan, we don't test... I suppose most labs score zero on The Joel Test). It is that we program, we program a lot and we still mostly use Fortran!

Now I'm not implying that Fortran is not a good language. Fortran is just fine. What I'm saying is that Fortran is a very bad choice for non-programmers who do a lot of programming.

Over the years some have turned to other languages like C, C++ and Matlab. Some, perhaps more enlightened, have even turned to Python. Naturally, I, too, turned to other languages and soon began wondering about what else might be out there (language-wise) to help a poor physicist do his job.

Could there be a language that won't break everything apart just because a number is consider "special"? (What's a float or a double? Why do I care?) Could there be one that implements a few-lines-of-equations-model in a few lines of code? One that will help me in my work not make the problem just ... different?

I've been googling and poking around on stackoverflow and programmers.SE for some answers but nothing particularly interesting came up. It, thus, seems that I have to try them all.

"Rules"
no extra packages
idiomatic code for each language
efficiency is not so much of an issue
extendable code


[TOC]

The system

I'll be solving, as an example, the following 2D system of ordinary differential equations:

\begin{align*} \frac{df_1(x)}{dx} =& 1 - f_2(x)\\ & \\
\frac{df_2(x)}{dx} =& f_1(x) - x
\end{align*}
with initial conditions:
\begin{align*}
x_0 &= 0\\
f_1(x_0) &= 1\\
f_2(x_0) &= 0
\end{align*}
in the range [0, 3]

The exact solution is f_1(x) = \cos{x} + x, f_2(x) = \sin{x}.

System solution

To numerically solve the system we treat the function values f_i(x) as parameters of the differential functions:
\begin{align*}
F_1(x,y_1,y_2) &= \frac{df_1(x,y_1,y_2)}{dx} =& 1 - y_2\\ && \\
F_2(x,y_1,y_2) &= \frac{df_2(x,y_1,y_2)}{dx} =& y_1 - x
\end{align*}

We will end up with a list of y_{i,0},\ldots, y_{i,n} values that will represent an approximation of the values f_i(x_k), x_k = x + k*\delta x.

The method: Runge Kutta-4

Runge Kutta method iterates from a current known value of f(x) to the "next" value at x+\delta x. So, knowing f(x_0) and choosing a small step \delta x the method gives us f(x_0 + \delta x) by calculating a "correction" to f(x_0).

We thus have:  f(x_0 + \delta x) = f(x_0) + \Delta() where the correction \Delta() depends on \delta x, x_0, f(x_0), f\prime(x,y) and is given by:

\begin{equation*}
\Delta(\delta x, x_0, f(x_0), f\prime(x,y)) = \frac{1}{6}\delta x (a + 2b + 2c + d)
\end{equation*}
where
\begin{align*}
a =& f\prime(x_0, f(x_0))\\
b = & f\prime(x_0 + \frac{1}{2}\delta x, f(x_0) + \frac{1}{2}a\delta x)\\
c = & f\prime(x_0 + \frac{1}{2}\delta x, f(x_0) + \frac{1}{2}b\delta x)\\
d = & f\prime(x_0 + \delta x, f(x_0) + c\delta x)\\
\end{align*}

We feed \Delta() with the initial known values.

In many dimensions

To apply RK4 to a system of ODEs

for each set of equations F_i = f_i\prime(x_0, y_1, \ldots, y_n) we have a different set of parameters a_i, b_i, c_i, d_i and a correction function f_i(x_0 + \delta x) = f_i(x_0) + \Delta(\delta x, a_i,b_i,c_i,d_i)

\begin{align*}
a_i &= F_i(x_0 ,& y_1,& \ldots ,& y_n) \\
b_i &= F_i(x_0+\frac{1}{2}\delta x ,& y_1 + \frac{1}{2}a_1\delta x ,& \ldots ,& y_n + \frac{1}{2}a_n\delta x)\\
c_i &= F_i(x_0+\frac{1}{2}\delta x ,& y_1 + \frac{1}{2}b_1\delta x ,& \ldots ,& y_n + \frac{1}{2}b_n\delta x)\\
d_i &= F_i(x_0+\delta x ,& y_1 + c_1\delta x ,& \ldots ,& y_n + c_n\delta x)\\
\end{align*}

 \Delta() = \frac{1}{6}\delta x (a_i + 2 b_i + 2 c_i + d)


f_i(x_0 + \delta x) = f_i(x_0) + \Delta(\delta x, a_i,b_i,c_i,d_i)

The Code

I'll use a C implementation as a reference.

C (reference)

Adapted from numerical recipes in C:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
void rk4(double y[], double dydx[], int n, double x, double h, double yout[], void (*derivs)(double, double [], double []))
{
    int i;
    double xh, hh, h6;
        double dym[n], dyt[n], yt[n];

    hh = h*0.5;
    h6 = h/6.0;
    xh = x+hh;

    for (i=0;i<n;i++) yt[i] = y[i] + hh*dydx[i];
    (*derivs)(xh,yt,dyt);

    for (i=0;i<n;i++) yt[i] = y[i] + hh*dyt[i];
    (*derivs)(xh,yt,dym);

    for (i=0;i<n;i++) {
        yt[i] = y[i] + h*dym[i];
        dym[i] = dym[i] + dyt[i];
    }

    (*derivs)(x+h,yt,dyt);

    for (i=0;i<n;i++)
        yout[i] =  y[i] + h6*(dydx[i] + dyt[i] + 2.*dym[i]);
}

void rkdump(double vstart[], int nvar, double t1, double t2, int nstep, void (*derivs)(double, double[], double[]))
{
    int i,k;
    double x,h;
    double v[nvar], vout[nvar], dv[nvar];

    for (i=0;i<nvar;i++) {
        v[i] = vstart[i];
    }
    x = t1;
    h = (t2-t1)/nstep;
    for (k=0;k<nstep;k++) {
        (*derivs)(x,v,dv);
        rk4(v,dv,nvar,x,h,vout,derivs);
        if ((double)(x+h) == x) printf("Step size too small in routine rkdump\n");
        x += h;
        for (i=0;i<nvar;i++) {
            v[i] = vout[i];
        }
        solution[0][k] = v[0];
        solution[1][k] = v[1];
    }
}

for example:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <complex.h>

#define N 2 //Dimensions
#define x0 0.0
#define x1 3.0
#define dt 0.01

double *solution[N] = {NULL};

void deriv1 (double x, double y[], double dydt[])
{
    dydt[0] = 1 - y[1];
    dydt[1] = y[0] - x;
}

int main (void)
{
    int i;
    int steps; // solution steps
    double starting_v[N]; // Starting vector

    steps = (x1-x0)/dt;

    for (i=0;i<N;i++)
        solution[i] = (double *)realloc((char *)solution[i], sizeof(double)*steps);

    /* -- initial values -- */
    starting_v[0] = 1.0;
    starting_v[1] = 0.0;

    rkdump(starting_v, N, x0, x1, steps, deriv1);

    for (i=0;i<steps;i++) {
        printf("%4f %4f %4f\n", (i+1)*dt, solution[0][i], solution[1][i]);
    }

    return 0;
}

Coming up: J, LISP, R, Fortran, Python, Haskell, Pari/GP, maxima.

Any other languages that you would like to see?

No tips yet.
Be the first to tip!
Like this post? Tip me with bitcoin!

1F17U6ucJtp4RcigS3RGgqtdS5b9MzPYDZ

Share on RedditShare on FacebookShare on Google+Tweet about this on TwitterPrint this page
{ Tags: , , , , }

1 comment

RSS / trackback

respond

  1. eelvex

    on July 29, 2011 at 23:07

    These great floating point tips are exactly the kind of things that you shouldn't have to know.

    [Reply]

Real Time Analytics