Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Maxima definition of Falkner-Skan governing equations
/* 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)) $
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.