/* 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))$