将 MATLAB 中的 ode45 代码转换为 python 代码

Question

如何使用 python 解决这个 MATLAB 颂歌问题这是 BC 的 IVP：

F'''+FF''-F'^2+1=0
F(0)=F'(0)=0,  F'(oo)=1

当前的matlab代码会生成如下图

和教科书上的解完全一样：

问题是我需要使用 python

重新编码相同的问题

% stagnation flow
clc; close all; clear all; clf;
tol=1e-3;
x=1; % f''(0)
dx=0.1;
Xf=3;
tspan=(0:dx:Xf);
Nt=Xf/dx+1;
for i=1:10000
  iter=i;
  x=x+0.0001;
  F = @(t,y)[-y(1)*y(3)-1+y(2)^2;y(1);y(2)];
  yo=[x;0;0];
  [t,y]= ode45(F,tspan,yo);
  y2=(y(Nt,2));
  % x=x/(y2^(3/2)) % f''())=f''(0)/[f'(inf)^(3/2)]
  if abs(y(Nt,2)-1.0)<tol, break, end
end
y2=(y(Nt,2));
% x=x/(y2^(3/2)) % f''())=f''(0)/[f'(inf)^(3/2)]

figure(1)
plot(t,y(:,1),t,y(:,2),t,y(:,3))
legend('fpp','fp','f');
xlabel('η=y(B/v)^2');
ylabel('f,fp,fpp');
title('Numerical solutions for stagnation flow')
fprintf ('η  \t\t     fp \n')
fprintf ('%.2f \t\t%.6f \n',t,y(:,2))

我曾尝试使用 python 对相同的问题进行编码，但我找不到任何关于此事的教程。

Answer 1

如果任务只是解决数学问题，可以说您已经在 Matlab 中“做错了”（在使用过于昂贵的方法的意义上）。当你想解决边界值问题时，你应该使用专用的 BVP 求解器 bvp4c，请参阅 Matlab 文档了解如何使用。

即使任务是实现单射法，寻根程序也应该升级为某种具有收敛阶数的方法，甚至二分法也应该比线性搜索更快。从 1, 1.1 开始使用未知二阶导数的正割法似乎也很有效。

在python中还有一个BVP解算器，如果可以的话，效果很好。

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

x0,xf = 0,3
F = lambda t,y: [-y[0]*y[2]-1+y[1]**2,y[0],y[1]]
bc = lambda y0,yf: [ y0[1], y0[2], yf[1]-1]

res = solve_bvp(F,bc,[x0,xf], [[0,0],[1,1],[0,xf-1]])
print(f"y''(0)={res.y[0,0]}")
x = np.linspace(x0,xf,150)
plt.plot(x,res.sol(x).T)
plt.plot(res.x,res.y.T,'x', ms=2)
plt.legend(["y''", "y'", "y"])
plt.grid(); plt.show()

导致阴谋

这个初始猜测仍然适用于 xf=20，但对 xf=30 无效，这可能需要更详细的初始猜测，初始曲线较短，线性部分较长。

对于更大的xf，以下初始化似乎效果很好

x = [0., 1.]
while x[-1] < xf: x.append(x[-1]*1.5)
xf = x[-1]
x = np.asarray(x); x[0]=0
y0 = x-x0-1; y0[0]=0
y1 = 0*x+1; y1[0]=0
y2 = 0*x; y2[0]=1
res = solve_bvp(F,bc,x, [y2,y1,y0], tol=1e-8)

Answer 2

试试这个可以让你的 matlab 代码运行在与 ODE45 相同的线性搜索逻辑下更快。我同意它太贵了。

clc; close all; clear all; clf;
flow = @(t, F) [-F(1)*F(3)-1+F(2)^2; F(1); F(2)];

tol = 1e-3;
x = 1;
y2 = 0; 

while abs(y2-1) >= tol
    [t, y] = ode45(flow, [0,3], [x;0;0]);
    x += 0.0001;
    y2 = y(end, 2);
end

plot(t, y)

这是 python

中的一个实现

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

def flow(t, F):
  return [-F[0]*F[2]-1+F[1]**2, F[0], F[1]]

tol = 1E-3
x = 1
y2 = 0

while np.abs(y2-1) >= tol:
  sol = solve_ivp(flow, [0,3], [x,0,0])
  x += 0.0001
  y2 = sol.y[1, -1]

plt.plot(sol.t, sol.y.T)
plt.legend([r"$F^{\prime \prime} \propto \tau $",r"$F^\prime \propto u$", r"$F \propto \Psi$"])
plt.axis([0, 3, 0, 2])
plt.xlabel(r'$\eta = y \sqrt{\frac{B}{v}}$')

这是一个响应与误差成正比的实现，运行s 更快。

clc; close all; clear all; clf;
flow = @(t, F) [-F(1)*F(3)-1+F(2)^2; F(1); F(2)];

tol = 1e-3;
x = 1;
error = 1; 

while abs(error) >= tol
    [t, y] = ode45(flow, [0,3], [x;0;0]);
    y2 = y(end, 2);
    error = y2 - 1;
    x -= 0.1*error;
end

plot(t, y)

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

def flow(t, F):
  return [-F[0]*F[2]-1+F[1]**2, F[0], F[1]]

tol = 1E-3
x = 1
error = 1

while np.abs(error) >= tol:
  sol = solve_ivp(flow, [0,3], [x,0,0])
  y2 = sol.y[1, -1]
  error = y2 - 1
  x -= 0.1 * error

plt.plot(sol.t, sol.y.T)
plt.legend([r"$F^{\prime \prime} \propto \tau $",r"$F^\prime \propto u$", r"$F \propto \Psi$"])
plt.axis([0, 3, 0, 2])
plt.xlabel(r'$\eta = y \sqrt{\frac{B}{v}}$')

将 MATLAB 中的 ode45 代码转换为 python 代码

Convert ode45 code from MATLAB to a python code

python

matlab

ode

ode45