# Numerically solving a partial differential equation in python with Runge Kutta 4

I'm supposed to solve the following partial differential equation in python using Runge-Kutta 4 method in time.

$$\frac{\partial}{\partial t}v(y,t)=Lv(t,y)$$

where $$L$$ is the following linear operator

$$L=\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)^{-1}\left\{-i\alpha(1-y^2)\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)-2i\alpha+\frac{1}{R}\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)^2\right\}$$

The spatial dimension is discretized along the $$N = 101$$ grid points from $$y = -1$$ to $$y = 1$$.

The boundary conditions are v[0] = 0 and v[N-1] = 0 and

$$\left.\frac{\partial v}{\partial y}\right|_{y=-1}=\left.\frac{\partial v}{\partial y}\right|_{y=1}=0$$

I have to discretize the second- and fourth-derivative operators using central finite difference and create a function for them. Using those functions, I can then create the $$L$$ operator.

Finally, using the initial condition $$v(y) = 0.02(1+\cos(\pi y))$$ time step $$dt= 0.01$$ and parameter values $$R=500$$ and $$\alpha=0.3$$, I have to find $$v$$ after 10s.

I have written the following code where the $$L$$ operator is put in a matrix form. The first few iterations seem to work giving somewhat reasonable results but then the values explode and I get an error. I'm not sure where the problem is, I might be misunderstanding how to implement the problem.

Another point I'm unsure about is whether the inverse of the $$\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)$$ operator can be understood as the inverse of the matrix representation of that operator or if that's not the case.

Can anybody give me some pointers as to how to solve this?

import numpy
import scipy
from scipy import linalg
from matplotlib import pyplot
%matplotlib inline

N=101
L=1.0

#creating the grid points
y=numpy.linspace(-L+0.04,L-0.04,N-4)
dy=(2*L)/(N-1)

v=numpy.empty(N-4)

# Operator of the second derivative acting on v with respect to y
def D2_v(N,dy):

# Setup the diagonal of the operator.
D = numpy.diag((-2.0/ (dy**2)) * numpy.ones(N-4))

# Setup the upper diagonal of the operator.
U = numpy.diag((1.0/ (dy**2)) * numpy.ones(N - 5), k=1)
# Setup the lower diagonal of the operator.
L = numpy.diag((1.0/ (dy**2)) * numpy.ones(N - 5), k=-1)
# Assemble the operator.
D2_v = D + U + L
# Setup the Neumann condition
D[-1, -1] = (-7.0/4.0) / (dy**2)
D[0, 0] = (-7.0/4.0) / (dy**2)

return D2_v

# Operator of the fourth derivative acting on v with respect to y
def D4_v(N,dy):

# Setup the diagonal of the operator.
D = numpy.diag((6.0/ (dy**4)) * numpy.ones(N-4))

# Setup the upper diagonal of the operator.
U = numpy.diag((-4.0/ (dy**4)) * numpy.ones(N - 5), k=1)
U2 = numpy.diag((1.0/ (dy**4)) * numpy.ones(N - 6), k=2)

# Setup the lower diagonal of the operator.
L = numpy.diag((-4.0/ (dy**4)) * numpy.ones(N - 5), k=-1)
L2 = numpy.diag((1.0/ (dy**4)) * numpy.ones(N - 6), k=-2)

# Assemble the operator.
D4 = D + U + L + U2 + L2

# Setup the Neumann condition
D4[-1, -1] = (5.0) / (dy**4)
D4[0, 0] = (5.0) / (dy**4)
D4[1, 0] = (-7.0/4.0) / (dy**4)
D4[-2, -1] = (-7.0/4.0) / (dy**4)

return D4

D2=D2_v(N,dy)
D4=D4_v(N,dy)

# Operator L
def L_v(N,y,dy,R,alpha):

Ma = alpha*numpy.identity(N-4)
My = numpy.diag(y)
I = numpy.identity(N-4)
Lv = numpy.linalg.inv(D2-Ma**2)*(-1j*Ma*(I-My**2)*(D2-Ma**2)+1j*
(-2.0)*Ma+(1.0/R)*(D4-D2*(Ma**2)-(Ma**2)*D2+(Ma**4)))

return Lv

#parameters
alpha=0.3
R=500.0
dt=0.01

L1=L_v(N,y,dy,R,alpha)

#initial value for v
for i in range(N-4):
v[i]=0.02*(1.0+numpy.cos(numpy.pi*y[i]))

#Runge Kutta implementation for 1000 time steps
for m in range(1000):

k1 = dt*numpy.matmul(L1,v)
k2 = dt*numpy.matmul(L1,v + 0.5*k1)
k3 = dt*numpy.matmul(L1,v + 0.5*k2)
k4 = dt*numpy.matmul(L1,v + k3)

v = v + (1.0/6.0)*(k1+2*k2+2*k3+k4)
• first question is how are you choosing the N. Are you sure that your scheme is stable with this choice? See CFL condition Dec 27 '19 at 5:58
• Try to use an implicit solver, semidiscretized problems are stiff and explicit method usually have really sever time-step restrictions.
– VoB
Dec 27 '19 at 8:27
• Hi, thanks for replying, this being for an assignment I have to use both the N = 101 value and method given above. I believe the problem lies in the implementation of the method or my understanding of how to solve the equation Dec 27 '19 at 12:51
• and thanks for the editing! Dec 27 '19 at 13:28
• I would use intermediate state vectors here for evaluation of the RHS of the system of ODEs, that would probably simplify the debugging; i.e., u={...}v, then solve (...)w=u, then d/dt(v) = w. Dec 27 '19 at 19:02