# Simple Runge-Kutta Scheme for nonlinear PDE

I am new to this community as well as to scientific programming. I programmed a simple 4th order Runge-Kutta for the 1-D Cahn-Hilliard Equation for some first simple calculations on pattern forming systems. It turned out to be extremely sensitive to the change of the size of the spatial grid. With 100 grid elements the time step needs to be like 100 times smaller in order to converge ... This is the code i used:

import numpy as np
import matplotlib.pyplot as plt

# Parameters

k_c = 1.
eps = 0.5
L=12. # Domainsize

size = 50 # Gridsize
dx = L/size

T = 100 # total time
dt = .001 # time step
n = int(T/dt)

# Random Initial Conditions

U = np.random.normal(0,0.01,size) # Normalverteilte Zufallswerte um 0

# Laplace-Operator

def laplacian(Z):
Zleft = np.roll(Z,1)
Zright = np.roll(Z,-1)
Zcenter = Z
return (Zright + Zleft -2*Zcenter)/dx**2

def DoubleLaplacian(Z):
Z2left = np.roll(Z,2)
Z2right = np.roll(Z,-2)
Zleft = np.roll(Z,1)
Zright = np.roll(Z,-1)
Zcenter = Z
return (Z2left - 4* Zleft + 6*Zcenter - 4*Zright + Z2right)/dx**4

# Equation Body

def equation(U):
return laplacian(U**3)-2*eps/(k_c**2)*laplacian(U)-eps*(k_c**4)*DoubleLaplacian(U)

# Runge-Kutta-Step

def RK4Step(U, dt, equation):
k1 = dt*equation(U)
U_temp = U + k1/2; k2 = dt*equation(U_temp);
U_temp = U + k2/2; k3 = dt*equation(U_temp);
U_temp = U + k3; k4 = dt*equation(U_temp);
U += (k1 + 2*k2 + 2*k3 + k4)/6

# Integration
for i in range(n):
RK4Step(U, dt, equation)


Since I am new to this kind of stuff I wanted to ask whether there is a conceptual fault in my thinking or coding or it is simply due to the nonlinearity. I would be as well pretty happy to get some hints how to improve the performance of the code.

• I have not revised your code, but googling the equation, I've found that it is a diffusion-like equation, and this probably answers for the scaling of the time-step with the thining of the grid. Have you tried to compute the CFL condition for this equation? – Hydro Guy Jul 8 '15 at 16:07
• i also don't have time for a code review at the moment, but about performance-if you're new to sci. programming it may be worth noting that you're using python, a high level language that manages memory for you. you'll typically have fewer lines of code, but you pay a price for ease of use. compared with a lower level code (e.g., c++) your python might be lots slower if you don't outsource computation heavy things to a lower level language and likely won't scale well. – aeroNotAuto Jul 8 '15 at 18:53

import cProfile, pstats, StringIO
pr = cProfile.Profile()
pr.enable()
# Integration
for i in range(n):
RK4Step(U, dt, equation)
pr.disable()
s = StringIO.StringIO()
sortby = 'cumulative'
ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
ps.print_stats()
print s.getvalue()


gives the following :

         27300018 function calls in 103.787 seconds

Ordered by: cumulative time

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
2    0.000    0.000  103.787   51.894 /Users/Leo/anaconda/lib/python2.7/site-packages/IPython/core/interactiveshell.py:3005(run_code)
1    0.466    0.466  103.787  103.787 <ipython-input-2-3cf9531a6a7d>:5(<module>)
100000    4.605    0.000  103.317    0.001 <ipython-input-1-3a5916441c4c>:44(RK4Step)
400000    9.648    0.000   98.712    0.000 <ipython-input-1-3a5916441c4c>:39(equation)
3200000   16.343    0.000   63.233    0.000 /Users/Leo/anaconda/lib/python2.7/site-packages/numpy/core/numeric.py:1279(roll)
800000   13.794    0.000   45.859    0.000 <ipython-input-1-3a5916441c4c>:23(laplacian)
400000   12.037    0.000   43.205    0.000 <ipython-input-1-3a5916441c4c>:29(DoubleLaplacian)
3200000   12.589    0.000   12.589    0.000 {numpy.core.multiarray.concatenate}
6400000   12.371    0.000   12.371    0.000 {numpy.core.multiarray.arange}


You can see how Laplacian and Double laplacian take about 88% of the computing time.

You can either rewrite your code in Fortran or C and wrap it with f2py or rewrite those functions in cython.

Another way is compiling your code using Nuitka, this decreases runtime to 64%.

You also see that roll is what takes most of the time. So you can start working on an optimised version of that.

         25600000 function calls in 64.852 seconds

Ordered by: cumulative time

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
3200000   17.494    0.000   64.852    0.000 /usr/local/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/core/numeric.py:1279(roll)
6400000   12.644    0.000   12.644    0.000 {numpy.core.multiarray.arange}
3200000   12.635    0.000   12.635    0.000 {numpy.core.multiarray.concatenate}
3200000    4.499    0.000   10.167    0.000 /usr/local/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/core/numeric.py:464(asanyarray)
3200000    6.472    0.000    6.472    0.000 {method 'take' of 'numpy.ndarray' objects}
3200000    5.668    0.000    5.668    0.000 {numpy.core.multiarray.array}
3200000    5.440    0.000    5.440    0.000 {method 'reshape' of 'numpy.ndarray' objects}


Finally, find a way to get your CFL number to define your $\Delta t$