python 刚性常微分方程系统的转换和求解

Question

>我有一个刚性的微分方程组给一阶ODE给出。这个系统是用枫树写的。Maple 使用的默认方法是 Rosenbrock 方法。现在我的任务是用python工具解决这些方程。

1)我不知道如何在python代码中编写方程式。

2)我不知道如何使用numpy，scipy，matplotlib或PyDSTool求解方程。对于库 PyDSTool，我根本没有找到任何示例，尽管我读到它非常适合僵硬系统。

法典：

import numpy
import scipy
import matplotlib
varepsilon = pow(10, -2); j = -2.5*pow(10, -2); e = 3.0; tau = 0.3; delta = 2.0
u0 = -math.sqrt(-1 + math.sqrt(varepsilon ** 2 + 12) / varepsilon) * math.sqrt(2) / 6
u = -math.sqrt(-1 + math.sqrt(varepsilon ** 2 + 12) / varepsilon) * math.sqrt(2) * (1 + delta) / 6
v = 1 / (1 - 2 / e) * math.sqrt(j ** 2 + (1 - 2 / e) * (e ** 2 * u ** 2 + 1))
y8 = lambda y1,y5,y7: 1 / (1 - 2 / y1) * math.sqrt(y5 ** 2 + (1 - 2 / y1) * (1 + y1 ** 2 * y7 ** 2))
E0 = lambda y1,y8: (1 - 2 / y1) * y8
Phi0 = lambda y1,y7: y1 ** 2 * y7

y08 = y8(y1=e, y5=j, y7=u0);
E = E0(y1=e, y8=y08); Phi = Phi0(y1=e, y7=u0)
# initial values
z01 = e; z03 = 0; z04 = 0; z05 = j; z07 = u0; z08 = y08;
p1 = -z1(x)*z5(x)/(z1(x)-2);
p3 = -z1(x)^2*z7(x);
p4 = z8(x)*(1-2/z1(x));
Q1 = -z5(x)^2/(z1(x)*(z1(x)-2))+(z8(x)^2/z1(x)^3-z7(x)^2)*(z1(x)-2);
Q3 = 2*z5(x)*z7(x)/z1(x);
Q4 = 2*z5(x)*z8(x)/(z1(x)*(z1(x)-2));
c1 = z1(x)*z7(x)*varepsilon;
c3 = -z1(x)*z5(x)*varepsilon;
C = z7(x)*varepsilon/z1(x)-z8(x)*(1-2/z1(x));
d1 = -z1(x)*z8(x)*varepsilon;
d3 = z1(x)*z5(x)*varepsilon;
B = z1(x)^2*z7(x)-z8(x)*varepsilon*(1-2/z1(x));
Omega = 1/(c1*d3*p3+c3*d1*p4-c3*d3*p1);
# differential equations
diff(z1(x), x) = z5(x);
diff(z3(x), x) = z7(x);
diff(z4(x), x) = z8(x);
diff(z5(x), x) = Omega*(-Q1*c1*d3*p3 - Q1*c3*d1*p4 + Q1*c3*d3*p1 + B*c3*p4 + C*d3*p3 + E*d3*p3 - Phi*c3*p4);
diff(z7(x), x) = -Omega*(Q3*c1*d3*p3 + Q3*c3*d1*p4 - Q3*c3*d3*p1 + B*c1*p4 - C*d1*p4 + C*d3*p1 - E*d1*p4 + E*d3*p1 - Phi*c1*p4);
diff(z8(x), x) = Omega*(-Q4*c1*d3*p3 - Q4*c3*d1*p4 + Q4*c3*d3*p1 + B*c1*p3 - B*c3*p1 - C*d1*p3 - E*d1*p3 - Phi*c1*p3 + Phi*c3*p1);
#features to be found and built curve
{z1(x), z3(x), z4(x), z5(x), z7(x), z8(x)}

Answer 1

上网漂流后，我原则上发现了一些东西：

import math
import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate
from scipy.signal import argrelextrema
from mpmath import mp, mpf
mp.dps = 50

varepsilon = pow(10, -2); j = 2.5*pow(10, -4); e = 3.0; tau = 0.5; delta = 2.0
u0 = -math.sqrt(-1 + math.sqrt(varepsilon ** 2 + 12) / varepsilon) * math.sqrt(2) / 6
u = -math.sqrt(-1 + math.sqrt(varepsilon ** 2 + 12) / varepsilon) * math.sqrt(2) * (1 + delta) / 6
v = 1 / (1 - 2 / e) * math.sqrt(j ** 2 + (1 - 2 / e) * (e ** 2 * u ** 2 + 1))
y8 = lambda y1,y5,y7: 1 / (1 - 2 / y1) * math.sqrt(y5 ** 2 + (1 - 2 / y1) * (1 + y1 ** 2 * y7 ** 2))
E0 = lambda y1,y8: (1 - 2 / y1) * y8
Phi0 = lambda y1,y7: y1 ** 2 * y7

y08 = y8(y1=e, y5=j, y7=u0);
E = E0(y1=e, y8=y08); Phi = Phi0(y1=e, y7=u0)
# initial values
z01 = e; z03 = 0.0; z04 = 0.0; z05 = j; z07 = u0; z08 = y08;
def model(x, z, varepsilon, E, Phi):
z1, z3, z4, z5, z7, z8 = z[0], z[1], z[2], z[3], z[4], z[5]
p1 = -z1*z5/(z1 - 2);
p3 = -pow(z1, 2) *z7;
p4 = z8*(1 - 2/z1);
Q1 = -pow(z5, 2)/(z1*(z1 - 2)) + (pow(z8, 2)/pow(z1, 3) - pow(z7, 2))*(z1 - 2);
Q3 = 2*z5*z7/z1;
Q4 = 2*z5*z8/(z1*(z1 - 2));
c1 = z1*z7*varepsilon;
c3 = -z1*z5*varepsilon;
C = z7*varepsilon/z1 - z8*(1 - 2/z1);
d1 = -z1*z8*varepsilon;
d3 = z1*z5*varepsilon;
B = pow(z1, 2)*z7 - z8*varepsilon*(1 - 2/z1);
Omega = 1/(c1*d3*p3+c3*d1*p4-c3*d3*p1);
# differential equations
dz1dx = z5;
dz3dx  = z7;
dz4dx  = z8;
dz5dx  = Omega*(-Q1*c1*d3*p3 - Q1*c3*d1*p4 + Q1*c3*d3*p1 + B*c3*p4 + C*d3*p3 + E*d3*p3 - Phi*c3*p4);
dz7dx  = -Omega*(Q3*c1*d3*p3 + Q3*c3*d1*p4 - Q3*c3*d3*p1 + B*c1*p4 - C*d1*p4 + C*d3*p1 - E*d1*p4 + E*d3*p1 - Phi*c1*p4);
dz8dx = Omega*(-Q4*c1*d3*p3 - Q4*c3*d1*p4 + Q4*c3*d3*p1 + B*c1*p3 - B*c3*p1 - C*d1*p3 - E*d1*p3 - Phi*c1*p3 + Phi*c3*p1);
dzdx = [dz1dx, dz3dx, dz4dx, dz5dx, dz7dx, dz8dx]
return dzdx
z0 = [z01, z03, z04, z05, z07, z08]
if __name__ == '__main__':
# Start by specifying the integrator:
# use ``vode`` with "backward differentiation formula"
r = integrate.ode(model).set_integrator('vode', method='bdf')
r.set_f_params(varepsilon, E, Phi)
# Set the time range
t_start = 0.0
t_final = 0.1
delta_t = 0.00001
# Number of time steps: 1 extra for initial condition
num_steps = np.floor((t_final - t_start)/delta_t) + 1
r.set_initial_value(z0, t_start)
t = np.zeros((int(num_steps), 1), dtype=np.float64)
Z = np.zeros((int(num_steps), 6,), dtype=np.float64)
t[0] = t_start
Z[0] = z0

k = 1
while r.successful() and k < num_steps:
r.integrate(r.t + delta_t)
# Store the results to plot later
t[k] = r.t
Z[k] = r.y
k += 1  
# All done!  Plot the trajectories:
Z1, Z3, Z4, Z5, Z7, Z8  = Z[:,0], Z[:,1] ,Z[:,2], Z[:,3], Z[:,4], Z[:,5]
plt.plot(t,Z1,'r-',label=r'$r(s)$')
plt.grid('on')
plt.ylabel(r'$r$')
plt.xlabel('proper time s')
plt.legend(loc='best')
plt.show()
plt.plot(t,Z5,'r-',label=r'$frac{dr}{ds}$')
plt.grid('on')
plt.ylabel(r'$frac{dr}{ds}$')
plt.xlabel('proper time s')
plt.legend(loc='best')
plt.show()

plt.plot(t, Z7, 'r-', label=r'$frac{dϕ}{ds}$')
plt.grid('on')
plt.xlabel('proper time s')
plt.ylabel(r'$frac{dϕ}{ds}$')
plt.legend(loc='upper center')
plt.show()

但是，回顾图书馆获得的解决方案， 我遇到了 scipy 和 Maple 获得的解决方案不一致的问题。问题的本质是解决方案快速振荡，枫树使用罗森布洛克的方法高精度地捕捉这些振荡。虽然 Pythonn 在使用向后微分方法时存在问题：

r = integrate.ode(model).set_integrator('vode', method='bdf')

http://www.scholarpedia.org/article/Backward_differentiation_formulas

我尝试了所有集成模式：">vode";"兹沃德";"苏打" ;"多普里5" ;"dop853"和我发现最适合的模式"vode"仍然不能满足我的需求...... https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html

因此，此方法捕获 j ~ 10^{-5}-10^{-3} 范围内的振荡。虽然枫树对任何j都显示出良好的结果。

我呈现 scipy 获得的结果 j ~ 10^{-2}：

在此处输入图像描述

在此处输入图像描述

以及枫树对j~10^{-2}得到的结果：

在此处输入图像描述

在此处输入图像描述

重要的是，振荡是物理解决方案！也就是说，Python 严重捕获了 j ~ 10^{-2}((.谁能告诉我我做错了什么？？如何看待集成的绝对误差？