odeint 没有给我预期的结果

Question

我从 ODE 开始，我有一个 python 代码，它首先尝试使用 odeint 函数求解 DOE，然后将该解与您自己计算的 ODE 的解进行比较（您必须将其插入代码中）。然而它们并不完全一致。

我试过软件程序，我手动计算的解决方案似乎没问题，所以我想知道为什么 odeint 函数没有给我预期的结果。

import numpy as np
import matplotlib.pyplot as plt

# Define a function which calculates the derivative
def dy_dx(y,x):
    dydx = 2.0*x*(y-1)/(1-x**2) #
    return dydx

xs = np.linspace(-4,4,100) 

for i in range(-10,10,1):
    y0 = i
    ys = odeint(dy_dx, y0, xs) 
    ys = np.array(ys).flatten() 
    plt.plot(xs, ys);

# SECOND PART:

y_exact =  1+(y0)/(1-x**2) # manually computed solution

y_difference = ys - y_exact 
plt.subplot(2, 1, 1)
plt.plot(xs, ys, xs, y_exact, "--"); 
plt.title("Difference between exact solution and computed solution");

所以我添加了 "range() part" 以查看它如何随着不同的初始条件而变化，但它们都沿着轴 x=-1

我手动找到的解决方案在那里有一个限制，但它不仅仅是一行，你可以看到如果你运行第二部分或访问https://www.symbolab.com/solver/ordinary-differential-equation-calculator/2x%5Cleft(y-1%5Cright)dx%2B%5Cleft(x%5E%7B2%7D-1%5Cright)dy%3D%200

我只是想知道错误在哪里，或者为什么 odeint 会给出这样的结果。difference between odeint and my solution

我还应该补充一点，ODE 可能有点奇怪，因为积分时会得到绝对值。这可能是相关的。

Answer 1

你的数值解的初始条件是y(-4)=y0，因为odeint取时间间隔的第一个点作为初始时间。相应地，您必须将确切的解决方案更改为

y_exact =  1+((y0-1)*(1-xs[0]**2))/(1-xs**2)

您可以使用您选择的工具进行检查，例如 WA y'(x)=2*x*(y(x)-1)/(1-x^2), y(a)=b。

正如您还观察到的那样，您的 ODE 在 x=-1 和 x=1 处是奇异的，因此任何解决方案都只有由这些点创建的 3 个子区间中的任何一个作为定义域，即一个包含初始时间。此外，数值方法不能很好地接近奇点，因此对于一些小但不太小的 delta，您必须将积分域限制为 [-4, -1-delta]。或者，如果您真的想探索初始时间为 0 的中间区间，则需要执行两次积分，一次从 0 到 -1+delta，一次从 0 到 1-delta.

如果考虑无实奇点的相关微分方程

y'(x) = -2*x*(y-1)/(1+x^2),  y(x0) = y0

有精确解

y(x) = 1 + (y0-1)*(1+x0^2)/(1+x^2)

通过修改后的代码实现

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define a function which calculates the derivative
def dy_dx(y,x): return  -2.0*x*(y-1)/(1+x**2) #

# Flow function
def phi(x,x0,y0): return 1 + (y0-1)*(1+x0**2)/(1+x**2)

xs = np.linspace(-4,4,100) 

for i in range(-10,10,2):
    y0 = 1+0.1*i
    ys = odeint(dy_dx, y0, xs) 
    ys = ys.flatten() 
    plt.plot(xs, phi(xs,xs[0],y0),c='lightgray', lw=4);
    plt.plot(xs, ys, lw=1, label='$y_0=%.2f$'%y0);
plt.grid(); plt.xlim(-5.5,4.2); plt.legend(); plt.show()

如果您发现（彩色）数值解很好地位于（底层灰色）精确解的中心位置，您会得到以下图。

odeint 没有给我预期的结果

odeint doesn't give me what expected

python-2.x

ode

odeint