Programming of 4th order Runge-Kutta in advection equation in python

Question

%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from math import pi

# wave speed
c = 1
# spatial domain
xmin = 0
xmax = 1
#time domain
m=500; # num of time steps 
tmin=0
T = tmin + np.arange(m+1);
tmax=500

n = 50 # num of grid points

# x grid of n points
X, dx = np.linspace(xmin, xmax, n+1, retstep=True);
X = X[:-1] # remove last point, as u(x=1,t)=u(x=0,t)

# for CFL of 0.1
CFL = 0.3
dt = CFL*dx/c


# initial conditions
def initial_u(x):
    return np.sin(2*pi*x)

# each value of the U array contains the solution for all x values at each timestep
U = np.zeros((m+1,n),dtype=float)
U[0] = u = initial_u(X);




def derivatives(t,u,c,dx):
    uvals = [] # u values for this time step
    for j in range(len(X)):
        if j == 0: # left boundary
            uvals.append((-c/(2*dx))*(u[j+1]-u[n-1]))
        elif j == n-1: # right boundary
            uvals.append((-c/(2*dx))*(u[0]-u[j-1]))
        else:
            uvals.append((-c/(2*dx))*(u[j+1]-u[j-1]))
    return np.asarray(uvals)


# solve for 500 time steps
for k in range(m):
    t = T[k];
    k1 = derivatives(t,u,c,dx)*dt;
    k2 = derivatives(t+0.5*dt,u+0.5*k1,c,dx)*dt;
    k3 = derivatives(t+0.5*dt,u+0.5*k2,c,dx)*dt;
    k4 = derivatives(t+dt,u+k3,c,dx)*dt;
    U[k+1] = u = u + (k1+2*k2+2*k3+k4)/6;

# plot solution
plt.style.use('dark_background')
fig = plt.figure()
ax1 = fig.add_subplot(1,1,1)

line, = ax1.plot(X,U[0],color='cyan')
ax1.grid(True)
ax1.set_ylim([-2,2])
ax1.set_xlim([0,1])
def animate(i):
    line.set_ydata(U[i])
    return line,

I want to program in Python an advection equation which is (∂u/∂t) +c (∂u/∂x) = 0. Time should be discretized with Runge-kutta 4th order. Spatial discretiziation is 2nd order finite difference. When I run my code, I get straight line which transforms into sine wave. But I gave as initial condition sine wave. Why does it start as straight line? And I want to have sine wave moving forward. Do you have any idea on how to get sine wave moving forward? I appreciate your help. Thanks in advance!

First, you need to use less global variables, esp. the `U` array. Pass more as parameters. Note that for RK4 you will need 3 derivative evaluations at states that are **not** Euler steps from the last time node. Also, the `k` vectors should be the result of the derivative function (this change might be enough to correct your code). To further simplify the code, think about how you could implement the second derivative using slices or convolution. — Lutz Lehmann, Jan 05 '20 at 21:58
Thanks for your answer. What do you mean exactly by 3 derivative evaluations at states that are not Euler step from the last time node? Can you clarify it? — alimuradpasa, Jan 09 '20 at 20:06
Reading again later, that was not relevant here. You are just using a very unconventional method to get the stage values. However, the point after that remains, the function `u` computes a point/state, while the `k` vectors should contain slopes/derivatives. — Lutz Lehmann, Jan 09 '20 at 20:18
I included Derivative function in k vectors. So k values are calculated by using Derivative function. But I get now strong oscilliations. I edited my code. So you can see its new form. — alimuradpasa, Jan 09 '20 at 20:28
You are doing something very strange in dividing by `(dx+deltax)`. By its use, `deltax` should rather be `deltat` and at no point it is added to `dx`. — Lutz Lehmann, Jan 09 '20 at 20:30
Ok. I understand. Where should I add deltat in Derivative function then? — alimuradpasa, Jan 09 '20 at 20:47
You don't. The derivative evaluation should not depend on dt or dt/2. — Lutz Lehmann, Jan 09 '20 at 20:51

Lutz Lehmann · Accepted Answer · 2020-01-10T11:49:39.653

0

While superficially your computation steps are related to the RK4 method, they deviate from the RK4 method and the correct space discretization too much to mention it all.

The traditional way to apply ODE integration methods is to have a function derivatives(t, state, params) and then apply that to compute the Euler step or the RK4 step. In your case it would be

def derivatives(t,u,c,dx):
    du = np.zeros(len(u));
    p = c/(2*dx);
    du[0] = p*(u[1]-u[-1]);
    du[1:-1] = p*(u[2:]-u[:-2]);
    du[-1] = p*(u[0]-u[-2]);
    return du;

Then you can do

X, dx = np.linspace(xmin, xmax, n+1, retstep=True);
X = X[:-1] # remove last point, as u(x=1,t)=u(x=0,t)

m=500; # number of time steps
T = tmin + np.arange(m+1);

U = np.zeros((m+1,n),dtype=float)
U[0] = u = initial_u(X);
for k in range(m):
    t = T[k];
    k1 = derivatives(t,u,c,dx)*dt;
    k2 = derivatives(t+0.5*dt,u+0.5*k1,c,dx)*dt;
    k3 = derivatives(t+0.5*dt,u+0.5*k2,c,dx)*dt;
    k4 = derivatives(t+dt,u+k3,c,dx)*dt;
    U[k+1] = u = u + (k1+2*k2+2*k3+k4)/6;

This uses dt as computed as the primary variable in the time stepping, then constructs the arithmetic sequence from tmin with step dt. Other ways are possible, but one has to make tmax and the number of time steps compatible.

The computation up to this point should now be successful and can be used in the animation. In my understanding, you do not produce a new plot in each frame, you only draw the graph once and after that just change the line data

# animate the time data
line, = ax1.plot(X,U[0],color='cyan')
ax1.grid(True)
ax1.set_ylim([-2,2])
ax1.set_xlim([0,1])

def animate(i):
    line.set_ydata(U[i])
    return line,

etc.

edited Jan 10 '20 at 11:49

answered Jan 09 '20 at 20:50

Lutz Lehmann

25,219
2
22
51

Thanks for your help! But I am getting error since t+0.5dt is float number. Error is "Only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices" which appears in k2 = derivatives(t+0.5*dt,u+0.5*k1,c,dx)*dt; – alimuradpasa Jan 09 '20 at 21:51
It is in general a bad idea to use names that normally imply a floating point variable as integer index. You need also to ensure that you use matching interfaces in the function calls, in my way of doing it the times are not actually used. – Lutz Lehmann Jan 09 '20 at 21:58
But in your way, you used t which is index . When 0.5dt is added to it, it will be floating point number as index. How should I solve this problem? Can you clarify please? – alimuradpasa Jan 09 '20 at 23:25
Thanks for your edited code. I get an error now in line k1 = derivatives(t,u,c,dx)*dt; error is like that - can't multiply sequence by non-int of type 'numpy.float64' – alimuradpasa Jan 10 '20 at 17:45
I'm not sure what you did, I got a nice smooth animation. If you use your derivative function or something similar, you need to convert the constructed list into an numpy array, `return np.asarray(uvals)`. – Lutz Lehmann Jan 10 '20 at 18:27
Yes, that should work. Note that the times are not correct, the time step in T is one, while in the simulation it is dt. Add in back the call to the function animation, and it should be complete. – Lutz Lehmann Jan 10 '20 at 22:57
I am getting still same error "IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices" in line k2 = derivatives(t+0.5*dt,u+0.5*k1,c,dx)*dt. – alimuradpasa Jan 13 '20 at 12:24
Replace the global `U[t-1]` with the local `u`, then you will also get no error message about non-integer `t`. – Lutz Lehmann Jan 13 '20 at 13:17
Thanks a lot. Yes it worked. But sine wave stays constant. It gives an warning which is RuntimeWarning: divide by zero encountered in double_scalars. How can I get sine wave moving forward? – alimuradpasa Jan 13 '20 at 16:28
Look how you are grouping the terms in the second derivative computation. You are dividing too much, it should only be by `2*dx`. – Lutz Lehmann Jan 13 '20 at 16:30
Yes. I corrected it. I get now no warning. Bu sine wave stays still constant. – alimuradpasa Jan 13 '20 at 16:43
Okay. I corrected animation function. So I now get nice sine wave. Thanks a lot! – alimuradpasa Jan 13 '20 at 16:51

Programming of 4th order Runge-Kutta in advection equation in python

1 Answers1