Julia programming

Revision as of 19:17, 10 March 2020

One can generate an array from a range 1:10 using collect():

collect(1:10)

10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

Let us now turn to a typical numerical problem, solving differential equations. We'll start with Ordinary Differential Equations (ODE). Namely, we consider equations of the type:

$$\tag{1} y'(t)=f(t,y(t))$$

The most famous algorithm is Euler's method and the most popular the RK4 variant of the Runge–Kutta methods. In between are a variety of possible implementations, many related in some way or the other. Methods can be explicit (the solution only depends on quantities already known) or implicit (the solution involves other unknowns which requires solving intermediate equations). The backward Euler method is an implicit method. Euler's method is the simplest one and the typical starting point to introduce the basic idea. Not surprisingly, it is also one of the least efficient. We discretize time into steps $t_1$, $t_2$, $\cdots$, $t_n$ separated by the same interval $h$. The solutions at these times are denoted $y_n=y(t_n)$. Back to $(1)$, approximating $y'(t)\approx (y(t+h)-y(t))/h$, we have:

$$y_{n+1}=y_n+hf(t_n,y_n)$$

This is very simple since from the solution $y_n$ at time $t_n$ we can compute the solution at $t_{n+1}$. Explicitly.

Let us solve, for instance, $y'=-y$, whose solution is the exponential function.

npts = 100; # number of points
y=0.0*collect(1:npts); # predefine the array
h=.1; # timesteps, so total time = npts*h=10
# The differential equation to solve
function f(t,y)
    -y
end
y[1]=1; # initial condition
# The numerical method (Euler's)
for i=1:npts-1
    y[i+1]=y[i]+h*f((i-1)*h,y[i])
end

And we can now compare the numerical solution to the real one:

plot(((1:npts).-1)*h, [y[i] for i=1:npts], label="Euler", xlabel="t", ylabel="y'=-y")
plot!(x->exp(-x),0, 10, label="Exact")

Let us try the more difficult example:

$$y'={y\over 2}+2\sin(3t)$$

whose solution can be found with an integrating factor [1]

$$y=-{24\over37}\cos(3t)-{4\over37}\sin(3t)+\left(y(0)+{24\over 37}\right)e^{t/2}$$

In this case the differential equation is defined as:

function f(t,y)
    (y/2)+2sin(3t)
end

The Volterra-Lotka model reads:

Euler's method in this case reads:

$$\begin{align} y_{1,n+1}&=y_{1,n}+hf_1(t_n,y_{1,n},y_{2,n})\\ y_{2,n+1}&=y_{2,n}+hf_2(t_n,y_{1,n},y_{2,n}) \end{align}$$

The functions are implemented as:

function f1(t,y1, y2)
    α*y1-β*y1*y2
end
function f2(t,y1, y2)
    δ*y1*y2-γ*y2
end

npts=10000
h=.001
 
y1=0.0*collect(1:npts);
y2=0.0*collect(1:npts);
 
y1[1]=1;
y2[1]=.1;
 
α=2/3;
β=4/3;
γ=1;
δ=1;
 
@time for i=1:npts-1
    y1[i+1]=y1[i]+h*f1((i-1)*h,y1[i],y2[i])
    y2[i+1]=y2[i]+h*f2((i-1)*h,y1[i],y2[i])
    end

plot([[y1[i] for i=1:npts], [y2[i] for i=1:npts]])

plot([(y1[i], y2[i]) for i=1:npts])

Before we jump to Runge-Kutta methods, let us look at fairly straightforwrad enhancement of Euler's method.

Backward Euler method, or Implicit Euler method,

Heun's Method is one of the simplest of a class of methods called predictor-corrector algorithms (Euler's method is only a predictor method).

We now cover one of the most powerful predictor-corrector algorithms of all, one which is so accurate that most computer algorithms use it as default: the fourth-order Runge-Kutta Method. RK methods are optimized predictor-corrector methods, so they push the idea of Heun's method to the top.