Convert ode45 code from MATLAB to a python code

Question

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

and it is identical to the textbook solution:

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.

Answer 1

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

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)

Answer 2

Try this to make your matlab code run faster under the same logic of linear search with ODE45. I agree it's too expensive.

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)

Here is an implementation in 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}}$')

And here is an implementation where the response is proportional to the error, which runs faster.

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}}$')