(* Install BVPh 1.0 for Mathematica 5.2. *) (* For Mathematica 8.0, please replace BVPh1_0.txt by BVPh1_1.txt *) < BVP in a finite domain [0,a] *) (* TypeEQ = 2 -> eigenvalue problem in a finite domain [0,a] *) (* TypeL = 1 -> Chebyshev polynomial as base function *) (* TypeL = 2 -> Hybrid-base approximation *) (* TypeBase = 1 -> sine + polynomial *) (* TypeBase = 2 -> cosine + polynomial *) (* ApproxQ = 0 -> do NOT approximate the function *) (* Appprox = 1 -> approximate the function *) (************************************************************************) TypeEQ = 1; ApproxQ = 0; ErrReq = 10^(-10); zRintegral = 10; H[z_] := z; (************************************************************************) (* Define the governing equation *) (************************************************************************) kappa = -1/3; alpha = 2/5; f[z_,u_,lambda_]:=alpha^2*D[u,{z,3}]+(1+kappa)/2*u*D[u,{z,2}]-kappa*D[u,z]^2; (***********************************************************************) (* Define Boundary conditions *) (***********************************************************************) zL = 1; zR = infinity; OrderEQ = 3; BC[1,z_,u_,lambda_] := Limit[u, z -> 1 ]; BC[2,z_,u_,lambda_] := Limit[D[u,z] - 1, z -> 1 ]; BC[3,z_,u_,lambda_] := Limit[D[u,{z,2}] - mu/alpha, z -> 1 ]; (************************************************************************) (* Define initial guess *) (************************************************************************) mu = 1; beta = (1+kappa)/(1-kappa); u[0] = Module[{temp}, temp[1] = (4-3*beta+mu/alpha)/(2-beta); temp[2] = (3-3*beta+mu/alpha)/(1-beta); temp[3] = (2-2*beta+mu/alpha)/(1-beta)/(2-beta); temp[1]*z^beta - temp[2]*z^(beta-1) + temp[3]*z^(2*beta-2) ]; (************************************************************************) (* Define output term *) (************************************************************************) output[z_,u_,k_]:= Print["output = ",alpha^2*D[u[k],{z,3}] /. z->1//N]; (************************************************************************) (* Defines the auxiliary linear operator *) (************************************************************************) L[u_] := Module[{temp}, temp[1] = 2*(2*beta-3); temp[2] = (beta-1)*(5*beta-6); temp[3] = 2*beta*(beta-1)^2; z^3*D[u,{z,3}] - temp[1]*z^2*D[u,{z,2}] + temp[2]*z*D[u,z] - temp[3]*u ]; (* Print input and control parameters *) PrintInput[u[z]]; Uexact = Module[{temp,t,t0,Ai,Bi,Ait,Bit}, Ai = AiryAi[t]; Bi = AiryBi[t]; Ait = D[Ai,t]; Bit = D[Bi,t]; Ait0 = Ait /. t -> t0; Bit0 = Bit /. t -> t0; temp[1] = Bit0*Ait - Ait0*Bit; temp[2] = Bit0*Ai - Ait0*Bi ; t0 = (Sqrt[6]*mu)^(-2/3); t = t0*(1 + mu*x); (36*mu)^(1/3)*temp[1]/temp[2]//Expand ]; Uzexact = D[Uexact,x]; Wz[k_] := Uz[k] /. z -> 1 + alpha*x; (* Set optimal c0 *) c0 = -8; Print[" c0 = ",c0]; (* Gain up to 10th-order HAM approximation *) BVPh[1,30];