pythonmatlabodeode45

Convert ode45 code from MATLAB to a python code

How can I solve this MATLAB ode problem using python This is the IVP with the BCs:

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

The current matlab code will generate the following plot enter image description here

and it is identical to the textbook solution: enter image description here

The problem is that I need to recode the same problem using 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))

I have tried to code the same problem using python but I couldn't find any tutorial about this matter.

Solution

  • If the task was just to solve the mathematical problem, one could say you already "did it wrong" in Matlab (in the sense of using a too expensive method). As you want to solve a boundary value problem, you should use the dedicated BVP solver bvp4c, see the Matlab documentation on how to.

    Even if the task was to implement a single-shooting approach, the root-finding procedure should be upgraded to some method with a convergence order, even bisection should be faster than the linear search. The secant method with the unknown second derivate starting at 1, 1.1 seems also to work well.

    In python there is also a BVP solver that works well if it works.

    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()

    resulting in the plot

    enter image description here

    This initial guess also still works with xf=20, but fails for xf=30, this probably needs a more detailed initial guess with a short initial curve and a long linear part.

    For larger xf the following intialization seems to work well

    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)