Блог

Как я могу решить эту проблему MATLAB ode с помощью python Это IVP с BCs:

 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 ('η  tt     fp n')
fprintf ('%.2f tt%.6f n',t,y(:,2))
 

Я пытался решить ту же проблему с помощью python, но не смог найти ни одного руководства по этому вопросу.

Если задача заключалась только в решении математической задачи, можно сказать, что вы уже «сделали это неправильно» в Matlab (в смысле использования слишком дорогого метода). Поскольку вы хотите решить граничную задачу, вам следует использовать специальный решатель BVP bvp4c , см. Документацию Matlab о том, как это сделать.

Даже если задача заключалась в реализации подхода с одним выстрелом, процедуру поиска корня следует обновить до некоторого метода с порядком сходимости, даже деление пополам должно быть быстрее, чем линейный поиск. Метод secant с неизвестным вторым производным, начинающимся с 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)
 

Попробуйте это, чтобы ваш код 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}}

И вот реализация, в которой ответ пропорционален ошибке, которая выполняется быстрее.

 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}}