blasiusSoluce.py
# -------------------------------------------------------------------------
#
# PYTHON for DUMMIES
# Probleme 8
#
# Solution detaillée
# Vincent Legat
#
# -------------------------------------------------------------------------
from numpy import *
# -------------------------------------------------------------------------
#
# -1- Equation de Blasius f''' + f f'' = 0
#
# dudt(t) = v(t)
# dvdt(t) = w(t)
# dwdt(t) = -u(t)*w(t)
#
#
def f(u):
dudt = u[1]
dvdt = u[2]
dwdt = - u[0] * u[2]
return array([dudt,dvdt,dwdt])
#
# -2- Méthode de bissection
#
def shoot(a,h,integrator):
X,U = integrator(0,[0,0,a],5,h,f)
return (1.0-U[-1,1])
def blasius(delta,nmax,tol,h,integrator):
a = delta[0]; fa = shoot(a,h,integrator)
b = delta[1]; fb = shoot(b,h,integrator)
n = 1; delta = (b-a)/2
if (fa*fb > 0) :
return 0,'Bad initial interval :-('
while (abs(delta) >= tol and n < nmax) :
delta = (b-a)/2; n = n + 1;
x = a + delta; fx = shoot(x,h,integrator)
print(" x = %14.7e (Estimated error %13.7e at iteration %d)" % (x,abs(delta),n))
if (fx*fa > 0) :
a = x; fa = fx
else :
b = x; fb = fx
if (n == nmax) :
return x,'Increase nmax : more iterations are needed :-('
return x,'Convergence observed :-)'