(* 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 = 2; TypeL = 2; TypeBase = 2; ApproxQ = 1; Ntruncated = 30; (************************************************************************) (* Define the governing equation *) (************************************************************************) L0[z_, u_] := Sqrt[1+z^2]*D[u,{z,2}] + Cos[Pi*z]*D[u,z]/z; F[z_,u_] := Exp[u]/(1+z^2) + (1+z)*Sin[u]; g[z_] := Sin[z^2 + Exp[-z]]; f[z_,u_,lambda_]:= L0[z,u] + lambda*F[z,u] - g[z]; (***********************************************************************) (* Define Boundary conditions *) (***********************************************************************) zR = Pi; OrderEQ = 2; BC[0,z_,u_,lambda_] := Limit[u - A, z->0] ; BC[1,z_,u_,lambda_] := Limit[D[u,z], z->0 ]; BC[2,z_,u_,lambda_] := u - D[u,z] - 3/5 /. z-> zR; A = 1/2; (************************************************************************) (* Define initial guess *) (************************************************************************) u[0] = (5*A+3)/10 + (5*A-3)/10*Cos[kappa*z]; kappa = 1; (************************************************************************) (* Define output term *) (************************************************************************) output[z_,u_,k_]:= Print["output = ", u[k] /. z -> zR //N ]; (************************************************************************) (* Defines the auxiliary linear operator *) (************************************************************************) omega[1] = Pi/zR; L[f_] := Module[{temp,numA,numB,i}, If[TypeL == 1, temp[1] = D[f,{z,OrderEQ}], numA = IntegerPart[OrderEQ/2]; numB = OrderEQ - 2* numA//Expand; temp[0] = D[f,{z,numB}]; For[i=1, i<=numA, i++, temp[1] = D[temp[0],{z,2}] + (kappa*omega[i])^2*temp[0]; temp[0] = temp[1]; ]; ]; temp[1]//Expand ]; (* Print input and control parameters *) PrintInput[u[z]]; (* Use 3rd-order HAM iteration approach *) iter[1,6,3];