After trying to use RK4 to solve the below system of equations, it appears the output had large and fast oscillations even with an adaptive time step I incorporated using the Cash-Karp method. I am now looking to existing solvers with more sophisticated tools (e.g. BDF methods) to solve PDEs, and came across odespy which I hope will work.
Here is the problem:
Tutorials for odespy are fairly detailed. Although helpful, an example of how to implement a system of PDEs with both time and spatial derivatives is missing, and I am wondering if odespy can actually solve such a problem. I am mainly curious since using the method of lines reduces this case to solving a system of a system of ODEs. Here is my code thus far:
import odespy
import matplotlib
matplotlib.use('pdf')
import os
import matplotlib.pyplot as plt
from pylab import *
import numpy as np
c = 3.*(10.**8.)
N = 500
zmax = 1000.
z = np.linspace(0.,zmax,N+1)
f_kwargs = dict(L=zmax, z=z)
y_1 = np.zeros((N+1),dtype=np.complex_) #should be [Nz x Nt] array
y_2 = np.zeros((N+1),dtype=np.complex_)
y_3 = np.zeros((N+1),dtype=np.complex_)
y_4 = np.zeros((N+1),dtype=np.complex_)
u = y_1,y_2,y_3,y_4
U_0 = np.zeros((4,N+1),dtype=np.complex_)
U_0[0][:] = np.cos(0.01)*(np.e**(z[:]/c))
U_0[1][:] = np.sin(0.01)*(np.e**(z[:]/c))
U_0[2][:] = 0.
U_0[3][:] = 0.
def rhs(u, t, L=None, beta=None, z=None):
N = len((u[0])[:]) - 1 #the PDE example specifies end condition, uses -1
dz = z[1] - z[0]
for i in range(0, N):
if i == 0:
OUT = [np.cos(0.01)*np.e**(-t),
np.sin(0.01)*np.e**(-t),
0.,
0.]
else:
OUT = [(y_2[i]*np.conj(y_3[i]) - y_1[i]),
(np.conj(y_3[i]*y_1[i]) - y_2[i]),
y_4[i],
(y_4[i] - ((y_3[i])**2.)*(1./dz)*np.conj(y_2[i+1] - y_2[i]))]
return OUT
solver = odespy.ForwardEuler(rhs, f_kwargs=f_kwargs,complex_valued=True)
solver.set_initial_condition(U_0)#solver.set_initial_condition(U_0)
dz = z[1] - z[0]
time_points = np.linspace(0.,600.,600) #t' = t/TR
dt = time_points[1] - time_points[0]
u, t = solver.solve(time_points)
This code does not work, and my main concerns are that I am:
1) Unable to properly formulate the equations in such a way odespy could handle a time derivative in my PDE (the 4th equation), and that the 3rd and 4th equations were derivatives in z as opposed to t (which is the derivative of the first two equations).
2) Unable to provide input to odespy such that it could take a system of PDEs, and then apply the method of lines approach outlined here. For example, one error already visible is my initial/boundary condition implementation. I used a 4 x (N+1) array for my output and initial/boundary conditions, but the solver is not expecting such an input from what I've gathered.
I'm still going through the documentation, but was wondering if anyone had any input for the above? Any suggestions at all on how to solve the above problem and comments on using odespy/numerical methods for a system of PDEs would be greatly appreciated.