(************************************************************************) (* Mathematica code for Blasius flow *) (* *) (* The optimal HAM is used to solve Blasius boundary-layer flow *) (* *) (* Dr. Shijun Liao *) (* Shanghai Jiao Tong University *) (* China *) (* June, 2010 *) (************************************************************************) < lambda*x ; fxx0[0] = D[f[0],{x,2}] /.x->0 ; lambda = N[4,100]; (************************************************************************) (* Define L *) (* This module defines the auxiliary linear operator *) (************************************************************************) L[f_] := D[f,{y,3}] + D[f,{y,2}]; (************************************************************************) (* Define Linv *) (* Find solution u of the following equation: L[u ] = f *) (************************************************************************) Linv[f_] := Module[{temp}, temp = DSolve[ L[w[y]] == f, w[y], y]; temp[[1,1,2]] /. C[_]->0 // Expand ]; (************************************************************************) (* the property of the inverse operator of L *) (* Linv[f_+g_] := Linv[f] + Linv[g] *) (************************************************************************) Linv[p_Plus] := Map[Linv,p]; Linv[c_*f_] := c*Linv[f] /; FreeQ[c,y]; (************************************************************************) (* Define delta[m] *) (* This module gives R[m] defined by (9.24) *) (************************************************************************) Getdelta[k_] := Module[{temp}, uyy[k] = D[u[k],{y,2}] //Expand ; uyyy[k] = D[uyy[k],y] //Expand ; temp[1] = Sum[u[k-n]*uyy[n],{n,0,k}]/2/lambda/lambda //Expand; temp[2] = uyyy[k] + temp[1] //Expand; delta[k] = Collect[temp[2], y^_.*Exp[_.]]; ]; (************************************************************************) (* Define GetErrorB[m] *) (* This module gives averaged residual error square *) (************************************************************************) GetErr[k_] := Module[{temp,dy,G,x,s}, Uyy[k] = D[U[k],{y,2}]; Uyyy[k] = D[Uyy[k],y]; error[k] = Uyyy[k] + U[k]*Uyy[k]/2/lambda/lambda; Ymax = 10; Nmax = 20; dy = N[Ymax/Nmax,100]; s = 0; For[j = 0, j <= Nmax, j = j + 1, x = j*dy; G[j] = error[k] /. y -> x //Expand; temp[2] = G[j]^2; s = s + temp[2]; ]; Err[k] = s/(Nmax+1) /. {AA->1-c1, BB-> 1-c2}; If[NumberQ[Err[k]],Err[k]//N//Print]; ]; (************************************************************************) (* Main Code *) (************************************************************************) ham[begin_,end_]:=Block[{uSpecial,B0,B2,temp,z,s}, time[0] = SessionTime[]; For[k=begin,k<=end,k=k+1, If[k==1, Print["-------------------------------------------------"] ]; If[k == 1 && OHAM == 0, Nstep = Infinity; If[NumberQ[c1], beta[1] = 1-c1, beta[1] = AA ]; If[NumberQ[c2], alpha[1] = 1-c2, alpha[1] = BB ]; Print["FINITE-parameter optimal HAM"]; Print[" c0 = ",c0]; Print[" c1 = ",c1]; Print[" c2 = ",c2]; If[c1 ==0 && c2 ==0, Print["This is the BASIC optimal HAM ! "] ]; ]; If[ k > 1 && OHAM == 0 , beta[k] = c1*beta[k-1]; alpha[k] = c2*alpha[k-1]; ]; If[k == 1 && OHAM > 0, If[Nstep == Infinity, Print["INFINITE-parameter optimal HAM"], Print["FINITE-parameter optimal HAM"] ]; alpha[1] = 1 - c2; c0 = 1; Print[" Nstep = ",Nstep]; Print[" c0 = ",c0]; Print[" c2 = ",c2]; ]; If[k > 1 && OHAM > 0, alpha[k] = c2*alpha[k-1] ]; If[k==1, Print["-------------------------------------------------"] ]; If[OHAM > 0 && k == 1, cc = {beta[1]} ]; If[OHAM > 0 && k > 1 && k <= Nstep, temp = Union[cc,{beta[k]}]; cc = temp ]; Print["k = ",k]; Getdelta[k-1]; temp = Linv[delta[k-1]]; S[k-1] = Collect[temp, y^_.*Exp[_.]]; temp = Sum[alpha[k-n]*u[n],{n,1,k-1}] ; uSpecial = temp + c0*Sum[beta[n]*S[k-n],{n,1,Min[k,Nstep]}]; B2 = D[uSpecial,y] /. y->0 ; B0 = -B2 - uSpecial /. y->0 ; temp = uSpecial + B0 + B2*Exp[-y] // Expand; u[k] = Collect[temp, y^_.*Exp[_.] ]; U[k] = Expand[U[k-1] + u[k]]; f[k] = U[k]/lambda /. {y -> lambda*x,AA->1-c1,BB->1-c2}; fx[k] = D[f[k],x]; fxx0[k] = D[fx[k],x] /. x->0 ; If[NumberQ[fxx0[k]], Print[" f''(0) = ", fxx0[k]//N ] ]; If[IntegerQ[k/5] && PRN == 1, time[k] = SessionTime[]; temp = time[k]-time[0]; Print["Used CPU times = ",temp, " (seconds) "]; ]; ]; Print["successful "]; ]; (************************************************************************) (* Define the parameters *) (************************************************************************) OHAM = 0; Nstep = Infinity; PRN = 1; c0 = -7/5; c1 = 0; c2 = 0; (* Gain the 15th-order homotopy-approximation *) ham[1,10];