(* 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; TypeL = 1; TypeBase = 2; ApproxQ = 0; (************************************************************************) (* Define the governing equation *) (************************************************************************) f[z_,u_,lambda_] := D[u,{z,4}] - beta * z * (1 + u^2); beta = 10 ; (***********************************************************************) (* Define Boundary conditions *) (***********************************************************************) zR = 1; OrderEQ = 4; alpha = 1/5 ; BC[1,z_,u_,lambda_] := Limit[u, z->0 ]; BC[2,z_,u_,lambda_] := D[u,z] /. z->zR; BC[3,z_,u_,lambda_] := D[u,{z,2}] /. z->zR; BC[4,z_,u_,lambda_] := Module[{temp}, temp[1] = D[u,{z,2}] /. z->0; temp[2] = D[u,{z,2}] /. z-> alpha; temp[1] - temp[2] //Expand ]; (************************************************************************) (* Define initial guess *) (************************************************************************) U[0] = u[0]; u[0] = sigma/(2*alpha-3)*((6*alpha-8)*z+6*(1-alpha)*z^2+2*alpha*z^3-z^4); sigma = .; (************************************************************************) (* Define output term *) (************************************************************************) output[z_,u_,k_]:= Print["output = ", N[ u[k] /. z->1, 24] ]; (************************************************************************) (* Defines the auxiliary linear operator *) (************************************************************************) 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]]; (* Gain 3rd-order HAM approximation *) BVPh[1,3]; (* Gain the squared residual of governing equation at 3rd-order approx. *) GetErr[3]; (* Gain optimal values of c0 and sigma *) res=Minimize[Err[3],{c0,sigma}]; Print["Minimum square residual = ", res]; (* Set c0 and sigma according to the above result *) c0=-1; sigma = 62/100; Print["c0 = ",c0," sigma = ",sigma]; (* Gain 10th-order HAM approximation *) BVPh[1,10];