Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Falkner-Skan IVP flat plate
import scipy as sp
from scipy.integrate import ode
from scipy.optimize import fsolve
from matplotlib import rc
#rc('text', usetex=True)
tick_size = 'large'
label_size = 'xx-large'
rc('xtick', labelsize=8)
rc('ytick', labelsize=8)
rc('legend', fontsize=12)
rc('axes', labelsize=12)
fig_width_pt = 469.75502
inches_per_pt = 1.0/72.27
golden_mean = (sp.sqrt(5.0)-1.0)/2.0
fig_width = fig_width_pt *inches_per_pt
fig_height = fig_width * golden_mean
fig_size = [fig_width, fig_height]
rc('figure', figsize=fig_size)
from matplotlib import pylab as plt
from fs import *
def int_fs(beta_0, beta, alpha, eta_inf):
args = sp.array([beta_0, beta, eta_inf])
y0 = sp.array([0.0, 0.0, alpha])
r = ode(fs.rates, fs.jac).set_integrator('zvode', rtol=1e-15, atol=1e-15, method='adams', with_jacobian=True, order=12)
t0 = 0.0
t1 = 1.0
nt = 100
dt = t1 / float(nt)
r.set_initial_value(y0, t0).set_f_params(args).set_jac_params(args)
for i in xrange(nt):
r.integrate(r.t+dt)
return(r)
def objective_func(x):
eta_inf = x[0]
alpha = x[1]
r = int_fs(1.0, 0.0, alpha, eta_inf)
return(sp.array([(r.y[1]-1.0).real, (r.y[2]).real]))
eta_inf = 1.0 # ~6
alpha = 1.0 # ~0.46960
x = sp.array([eta_inf, alpha])
foo = fsolve(objective_func, x, xtol=1e-15)
nt = 100
args = sp.array([1.0, 0.0, foo[0]])
y0 = sp.array([0.0, 0.0, foo[1]])
r = ode(fs.rates, fs.jac).set_integrator('zvode', rtol=1e-15, atol=1e-15, method='adams', with_jacobian=True, order=12)
t0 = 0.0
t1 = 1.0
nt = 100
dt = t1 / float(nt)
t = sp.linspace(0,1.0,nt+1)
y = sp.zeros((3,nt+1), dtype=float)
y[:,0] = y0.copy()
r.set_initial_value(y0, t0).set_f_params(args).set_jac_params(args)
for i in xrange(nt):
r.integrate(r.t+dt)
t[i+1] = r.t
y[:,i+1] = r.y.copy().real
plt.figure()
plt.plot(foo[0]*t,y[0,:])
plt.plot(foo[0]*t,y[1,:])
plt.plot(foo[0]*t,y[2,:])
plt.axis([0,foo[0],0,y.max()])
plt.xlabel('\eta')
plt.ylabel("f,f',f''")
plt.savefig("fs.png")
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.