Maxima definition of Falkner-Skan governing equations

falknerskan.mac

/* Falkner-Skan ODE

    Zhang, J., Chen, B., 'An iterative method for solving the
    Falkner-Skan equation'

    Nachtsheim, P., Swigert, P., 'Satisfaction of Asymptotic Boundary
    Conditions in Numerical Solution of Systems of Nonlinear Equations
    of Boundary-Layer Type,' NASA TN D-3004, Lewis Research Center,
    Cleveland, Oct, 1965.

*/

depends([f,u,v],xi)$

FS : diff(f,xi,3)/eta**3 + b_0*f*diff(f,xi,2)/eta**2 + b*(1-(diff(f,xi)/eta)**2) = 0 $

bc1_xi_1 : diff(f,xi) = eta $
bc2_xi_1 : diff(f,xi,2) = 0 $
bc3_xi_0 : diff(f,xi,2) = eta**2 * alpha $

sub_1 : diff(f,xi) = u*eta $
sub_2 : diff(f,xi,2) = v*eta**2 $
sub_3 : psubst([sub_1,sub_2], solve(FS,diff(f,xi,3))[1]) $
unknowns : [f,u,v] $
rates : diff(unknowns, xi) $
eqns : [sub_1,
        solve(subst(sub_2, diff(sub_1,xi)), rates[2])[1],
        solve(subst(sub_3, diff(sub_2,xi)),rates[3])[1]] $
A_aug : augcoefmatrix(eqns, unknowns) $

jac[i,j] := diff(rhs(eqns[i]), unknowns[j]) $
J : genmatrix(jac,3,3) $ 
eig_J : eigenvalues(J) $

/* dump some tex for documentation */
load("mactex-utilities") $
set_tex_environment_default("", "") $
texput(xi, "\\xi") $
texput(b, "\\beta") $
texput(b_0, "\\beta_0") $
texput(eta, "\\eta_{\infty}") $
fname : "J.tex" $
file_output_append : false $
with_stdout(fname, tex(J)) $
fname : "eqn_1.tex" $
with_stdout(fname, tex(eqn_1)) $
fname : "eqn_2.tex" $
with_stdout(fname, tex(eqn_2)) $
fname : "eqn_3.tex" $
with_stdout(fname, tex(eqn_3)) $
fname : "bc_1.tex" $
with_stdout(fname, tex(bc1_xi_1)) $
fname : "bc_2.tex" $
with_stdout(fname, tex(bc2_xi_1)) $
fname : "bc_3.tex" $
with_stdout(fname, tex(bc3_xi_0)) $

/* make substitutions so it's easy to use standard integrators */
int_subs : [f = y[1], u = y[2], v = y[3], b_0=arg[1], b=arg[2], eta=arg[3]] $

/* dump fortran for number crunching */
load(f90) $
fname : "f.f90" $
file_output_append : false $
F : transpose(psubst(int_subs, map(rhs,eqns))) $
with_stdout(fname, f90(F)) $
fname : "J.f90" $
file_output_append : false $
pJ : psubst(int_subs, J) $
with_stdout(fname, f90(pJ)) $