(* 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}] + u*D[u,{z,2}]/2; temp[2] = D[u,t]*D[u,{z,2}] - D[u,z]*D[u,z,t]; temp[1] + (1-t)* temp[2] // Expand ]; (***********************************************************************) (* Define Boundary conditions *) (***********************************************************************) zR = infinity; OrderEQ = 3; BC[1,z_,u_,lambda_] := u /. z-> 0 ; BC[2,z_,u_,lambda_] := D[u,z] - t /. z -> 0 ; BC[3,z_,u_,lambda_] := Limit[D[u,z], z -> zR ]; (************************************************************************) (* Define initial guess *) (************************************************************************) u[0] = t*(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->1}//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; ut[k] = D[u[k],t]//Expand; uzt[k] = D[uz[k],t]//Expand; uzzu[k] = Sum[uzz[i]*u[k-i],{i,0,k}]//Expand; uzuzt[k] = Sum[uz[i]*uzt[k-i],{i,0,k}]//Expand; uzzut[k] = Sum[uzz[i]*ut[k-i],{i,0,k}]//Expand; temp[1] = uzzz[k] + uzzu[k]/2 //Expand; temp[2] = (1-t)*(uzzut[k] - uzuzt[k]); delta[k] = temp[1] + temp[2]//Expand; ]; (* Print input and control parameters *) PrintInput[u[z,t]]; (* Set c0 *) c0 = -1/2; Print["c0 = ",c0]; (* gain the 10th-order HAM approximation by means of c0 =-1/2 *) 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]; ];