2

我正在尝试求解三阶非线性微分方程。我试图改变它，我得到了这个问题，这是一个二阶问题：

要解决的主要问题

我正在尝试实现四阶 Range-Kutta 算法，以便通过这样编写来解决它：

龙格-库塔问题

这是我的 Range-Kutta 算法代码：

import numpy as np
import matplotlib.pyplot as plt

''''X,Y = integrate(F,x,y,xStop,h).
4th-order Runge-Kutta method for solving the initial value problem {y}' = {F(x,{y})}, where {y} = {y[0],y[1],...,y[n-1]}.
x,y = initial conditions
xStop = terminal value of x 
h = increment of x used in integration
F = user-supplied function that returns the 
array F(x,y) = {y'[0],y'[1],...,y'[n-1]}.
'''

def integrate(F,x,y,xStop,h):
    
    def run_kut4(F,x,y,h):
        K0 = h*F(x,y)
        K1 = h*F(x + h/2.0, y + K0/2.0)
        K2 = h*F(x + h/2.0, y + K1/2.0)
        K3 = h*F(x + h, y + K2)
        return (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
    
    X =[]
    Y =[]
    X.append(x)
    Y.append(y)
    while x < xStop:
        h = min(h,xStop - x)
        y = y + run_kut4(F,x,y,h)
        x = x + h
        X.append(x)
        Y.append(y)
    return np.array(X),np.array(Y)

它适用于其他微分方程。

在这种情况下，函数 F 定义为：

F函数

主要代码是：

def F(x,y):
    F = np.zeros(2)
    F[0] = y[1]
    F[1] = (2*(1-x)/x**3)*y[0]**(-1/2)
    return F

x = 1.0
xStop = 20
y = np.array([0,0])
h = 0.2
X,Y = integrate(F,x,y,xStop,h)
plt.plot(X,Y)
plt.grid()
plt.show()

不幸的是，我收到了这个错误：

<ipython-input-8-8216949e6888>:4: RuntimeWarning: divide by zero encountered in power
  F[1] = (2*(1-x)/x**3)*y[0]**(-1/2)
<ipython-input-8-8216949e6888>:4: RuntimeWarning: divide by zero encountered in double_scalars
  F[1] = (2*(1-x)/x**3)*y[0]**(-1/2)

这与函数的初始值为0的事实有关，但我不知道如何摆脱它以便再次简化问题......

有人可以帮我找到其他选择吗？

谢谢您的帮助，

4

2 回答 2

1

你y是[0,0]并且在分母中y[0]**(-1/2)有除法运算，0这会给出 ZeroDivision 警告，并且在 double_scalars 中遇到的无效值是由于表达式y[0]**(-1/2)更改为NaN. 但是，这些是警告并且F正在返回 value array([ 0., nan])。您需要替换y[0]**(-1/2)，因为零的负幂是未定义的，或者如果适合您的需要，您可以使用接近零的极小值。也许您的方程在 (1,0) 处不连续。

于 2020-12-14T13:44:20.437 回答
0

您可以测试接近初始点的近似解，这些解是(1-x)该差的幂或幂级数。在最简单的情况下，您会y(x)=(1-x)^2得到x <= 1. 要获得解决方案，x>1您需要在平方根中取另一个符号。总的来说，这给出了远场行为x(t)=1+c*exp(-t)。

现在你可以遵循两种策略，

  • 将某个初始点的简化方程x=1-h与y(1-h)=h^2,y'(1-h)=-2h以及任意时间积分一起积分dt/dx=y(x)^(-1/2)，t(1-h)或
  • 将某个时间的原始方程t=T与 的导数之后的条件进行积分x(t)=1+c*exp(T-t)，即x(T)=1+c任意小c的 , x'(t)=-c, x''(T)=c. 理论上它们应该是相同的，实际上与固定步长RK4在两种情况下都会有差异。

切入所有方法。接近渐近线x(t)=1的解是单调的，因此时间可以用 表示x-1。这意味着导数可以（可能）表示为幂级数x-1

x'  (t) = c_1 * (x-1) + c_2 * (x-1)^2 + ...
x'' (t) = c_1 * x'(t) + 2c_2 * (x-1)*x' + ...
        = (c_1 + 2c_2*(x-1)+...)*(c_1+c_2*(x-1)+..)*(x-1)
        = c_1^2*(x-1)+3c_1c_2*(x-1)^2 + ...
x'''(t) = (c_1^2+6c_1c_2*(x-1)+...)*(c_1+c_2*(x-1)+..)*(x-1)
        = c_1^3*(x-1) + 7c_1^2c_2*(x-1)^2 + ...

and 

(1-x)/x^3 = -(x-1)*(1+(x-1))^(-3)=-(x-1)+3*(x-1)^2 + ...

so equating coefficients

c_1^3=-1 ==> c_1 = -1
7c_1^2c_2 = 3 ==> c_2 = 3/7

Given x(T) close enough to 1, the other initial values have to be 
x'(T)=-(x(T)-1) + 3/7*(x(T)-1)^2
x''(T)=x(T)-1 -9/7*(x(T)-1)^2

那么由于前两项的远场近似是

x'(t) = -(x-1) + 3/7 * (x-1)^2

Substitute u(t) = (x-1)^(-1) - 3/7

u'(t) = u(t)

(x(t)-1)^(-1) - 3/7 = ((x(T)-1)^(-1) - 3/7) * exp(t-T)

x(t) = 1 + (x(T)-1)*exp(T-t) / ( 1 - 3/7*(x(T)-1)*(1-exp(T-t)) )
于 2020-12-15T09:53:40.077 回答