计算科学 - 为 PDE 系统实现 odespy - 吾爱随笔录

为 PDE 系统实现 odespy

计算科学 pde Python 颂刚性耦合

2021-12-25 12:05:19

在尝试使用 RK4 求解以下方程组后，即使我使用 Cash-Karp 方法结合了自适应时间步长，输出似乎也具有大而快的振荡。我现在正在寻找具有更复杂工具（例如 BDF 方法）的现有求解器来求解 PDE，并且遇到了我希望能够工作的odespy 。

这是问题所在：

odespy 的教程相当详细。虽然有帮助，但缺少如何实现具有时间和空间导数的 PDE 系统的示例，我想知道 odespy 是否真的可以解决这样的问题。我主要是好奇，因为使用线的方法将这种情况减少到解决 ODE 系统的系统。到目前为止，这是我的代码：

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)

此代码不起作用，我主要担心的是：

1) 无法以 odespy 可以处理我的 PDE（第 4 个方程）中的时间导数的方式正确制定方程，并且第 3 和第 4 个方程是 z 中的导数，而不是 t（这是第一个的导数）两个方程）。

2) 无法向 odespy 提供输入，使其可以采用 PDE 系统，然后应用此处概述的线方法。例如，一个已经可见的错误是我的初始/边界条件实现。我为我的输出和初始/边界条件使用了一个 4 x (N+1) 数组，但求解器并不期望从我收集的内容中得到这样的输入。

我仍在阅读文档，但想知道是否有人对上述内容有任何意见？任何关于如何解决上述问题的建议以及对 PDE 系统使用 odespy/数值方法的评论都将不胜感激。

1个回答

Odespy需要一阶 ODE 系统，即形式为

\begin{aligned} (1) & \frac{d y}{d t} = f (y, t), \end{aligned}

$\begin{align} \frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y},t), \tag{1} \end{align}$

其中，对于您的特定问题，是仅依赖于时间的单变量复值函数（这意味着您需要首先通过离散空间维度来达到半离散形式，这就是方法的含义行），并且。请注意，定义的占位符变量（用于原始时空问题）。 $y_i$ $i=1 \ldots3$ $y_4$ $y_4 := \partial y_3 / \partial z$

所以如果我们在你原来的公式中去掉，我们得到时空系统： $y_4$

\begin{array}{rcl} \frac{\partial y_{1}}{\partial t} & = & (as before) \\ \frac{\partial y_{2}}{\partial t} & = & (as before) \\ \frac{\partial y_{3}}{\partial t} & = & \frac{\partial^{2} y_{3}}{\partial z^{2}} - \frac{\partial y_{3}}{\partial z} + y_{3}^{2} \frac{\partial y_{2}^{*}}{\partial z} . \end{array}

$\begin{eqnarray} \frac{\partial y_1}{\partial t} &=& \mbox{(as before)} \\ \frac{\partial y_2}{\partial t} &=& \mbox{(as before)} \\ \frac{\partial y_3}{\partial t} &=& \frac{\partial^2 y_3}{\partial z^2} - \frac{\partial y_3}{\partial z} + y_3^2 \frac{\partial y_2^*}{\partial z}. \end{eqnarray}$

所以剩下要做的是：

使用合适的有限差分公式对项)中的 RHS 进行离散化，以获得中的二阶和一阶偏导数，以得到类似于上述 (1) 的半离散系统，以及 $z$ $z$
将半离散系统传递给 Odespy（确保边界条件正确）。

确保选择可以处理系统时间尺度和其他属性的时间积分方法（例如，正向欧拉对于刚性系统来说效果不佳）。

后来编辑：仔细观察了系统，看起来是实变量的实值函数（初始/边界条件和右侧都是实数）。因此，如果您知道您的解决方案具有非零虚部，您应该仔细检查您的 PDE 系统。 $y_i$

其它你可能感兴趣的问题

上一篇选择具有已知导数的廉价黑盒函数全局优化的求解器/软件下一篇求解未知常数的积分方程