计算科学 - 使用 Runge Kutta 4 在 python 中数值求解偏微分方程 - 吾爱随笔录

我应该使用Runge-Kutta 4方法及时解决python中的以下偏微分方程。

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

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

在哪里 $L$ 是以下线性算子

L = {(\frac{\partial^{2}}{\partial y^{2}} - α^{2})}^{- 1} {- i α (1 - y^{2}) (\frac{\partial^{2}}{\partial y^{2}} - α^{2}) - 2 i α + \frac{1}{R} {(\frac{\partial^{2}}{\partial y^{2}} - α^{2})}^{2}}

$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\}$

空间维度沿 $N = 101$ 网格点从 $y = -1$ 到 $y = 1$ .

边界条件是v[0] = 0和v[N-1] = 0和

{\frac{\partial v}{\partial y} |}_{y = - 1} = {\frac{\partial v}{\partial y} |}_{y = 1} = 0

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

我必须使用中心有限差分对二阶和四阶导数运算符进行离散化，并为它们创建一个函数。使用这些函数，我可以创建 $L$ 操作员。

最后，使用初始条件 $v(y) = 0.02(1+\cos(\pi y))$ 时间步长 $dt= 0.01$ 和参数值 $R=500$ 和 $\alpha=0.3$ , 我必须找到 $v$ 10 秒后。

我编写了以下代码 $L$ 运算符以矩阵形式放置。前几次迭代似乎可以给出一些合理的结果，但随后值爆炸，我得到一个错误。我不确定问题出在哪里，我可能误解了如何解决问题。

我不确定的另一点是 $\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)$ 运算符可以理解为该运算符的矩阵表示的逆，或者如果不是这种情况。

任何人都可以给我一些关于如何解决这个问题的指示吗？

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)