(* 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); NgetErr = 100; zRintegral = 10; (************************************************************************) (* Define the governing equation *) (************************************************************************) f[z_,u_,lambda_] := Module[{temp}, temp[1] = D[u,{z,3}] + (1-t)*z/2*D[u,{z,2}] - t*(1-t)*D[u,z,t]; temp[2] = u*D[u,{z,2}] - D[u,z]^2; temp[1] + t* temp[2] // Expand ]; (***********************************************************************) (* Define Boundary conditions *) (***********************************************************************) zR = infinity; OrderEQ = 3; BC[1,z_,u_,lambda_] := Limit[u, z -> 0]; BC[2,z_,u_,lambda_] := Limit[D[u,z] - 1, z -> 0]; BC[3,z_,u_,lambda_] := Limit[D[u,z] , z -> zR ]; (************************************************************************) (* Define initial guess *) (************************************************************************) u[0] = 1 - Exp[-z]; (************************************************************************) (* Defines the auxiliary linear operator *) (************************************************************************) L[u_] := D[u,{z,3}] - D[u,z]; (************************************************************************) (* Define output term *) (************************************************************************) output[z_,u_,k_]:= Print["output = ",D[u[k],{z,2}]/.{z->0,t->0}//N]; (* Define Getdelta[k] *) Getdelta[k_]:=Module[{temp,i}, uz[k] = D[u[k],z]//Expand; uzz[k] = D[uz[k],z]//Expand; uzzz[k] = D[uzz[k],z]//Expand; uzuz[k] = Sum[uz[i]*uz[k-i],{i,0,k}]//Expand; uzzu[k] = Sum[uzz[i]*u[k-i],{i,0,k}]//Expand; uzt[k] = D[uz[k],t]//Expand; temp[1] = uzzz[k] + (1-t)*z/2*uzz[k]-t*(1-t)*uzt[k]//Expand; temp[2] = t*(uzzu[k] - uzuz[k]); delta[k] = temp[1] + temp[2]//Expand; ]; (* Print input and control parameters *) PrintInput[u[z,t]]; (* Set convergence-control parameter c0 *) c0 = -1/4; Print["c0 = ",c0]; (* Gain 6th-order HAM approximation *) Print[" c0 = ",c0]; BVPh[1,10]; (* Calculate the squared residual *) For[k=2, k<=10, k=k+2, GetErr[k]; err[k] = Integrate[Err[k],{t,0,1}]; Print[" k = ", k, " Squared residual = ", err[k]//N]; ];