(************************************************************************) (* BVPh 1.0 *) (* A Mathematica package for nonlinear boundary-value problems *) (* governed by a nonlinear ODE or some nonlinear PDEs. *) (* *) (* Version 1.0 is for Mathematica 5.2 *) (* *) (* Dr. Shijun Liao *) (* Shanghai Jiao Tong University *) (* China *) (* August, 2010 *) (************************************************************************) <0 // TrigReduce ]; (************************************************************************) (* the property of the inverse operator of L *) (* Linv[f_+g_] := Linv[f] + Linv[g] *) (************************************************************************) Linv[p_Plus] := Map[Linv,p]; Linv[c_*f_] := c*Linv[f] /; FreeQ[c,z]; (************************************************************************) (* Define Getdelta[k] *) (* Get delta[k] automatically based on the definition of f[z,u,lambda] *) (************************************************************************) Getdelta[k_]:=Module[{temp,phi,LAMBDA,w,lamb,q,eq,m,n,coeff,Coeff}, eq = f[z,phi[z,q],LAMBDA[q]]; temp[0] = D[eq,{q,k}]/k! // Expand; temp[1] = temp[0]/.D[phi[z,q],{z,m_},{q,n_}]->n!*diff[w[n],{z,m}]; temp[2] = temp[1]/.D[LAMBDA[q],{q,n_}]->n!*lamb[n]; temp[3] = temp[2]/.{phi[z,q]->w[0],LAMBDA[q]->lamb[0]}; temp[4] = temp[3]/.diff[w[m_],{z,0}]->w[m]; temp[5] = temp[4]/.{w->u,diff->D}; delta0[k] = temp[5] /. lamb->lambda; If[ApproxQ == 0, temp[-1] = delta0[k]//Expand//TrigReduce ]; If[ApproxQ == 1 && TypeL == 1, If[TypeEQ == 1, temp[-1] = ChebyApproxA[delta0[k],z,0,zR,Ntruncated], For[m=0,m<=k,m++,coeff[m]=Coefficient[temp[5],lamb[m]]]; coeff[-1] = temp[5]/.lamb[_]->0; For[n=-1,n<=k,n++, Coeff[n] = ChebyApproxA[coeff[n],z,0,zR,Ntruncated]; ]; lamb[-1] = 1; temp[-1] = Sum[Coeff[n]*lamb[n],{n,-1,k}]/.lamb->lambda; ]; ]; If[ApproxQ == 1 && TypeL == 2, If[TypeEQ == 1, temp[-1] = TrigApprox[delta0[k],z,zR,Ntruncated,TypeBase], For[m=0,m<=k,m++,coeff[m]=Coefficient[temp[5],lamb[m]]]; coeff[-1] = temp[5]/.lamb[_]->0; For[n=-1,n<=k,n++, Coeff[n] = TrigApprox[coeff[n],z,zR,Ntruncated,TypeBase]; ]; lamb[-1] = 1; temp[-1] = Sum[Coeff[n]*lamb[n],{n,-1,k}]/.lamb->lambda; ]; ]; delta[k] = temp[-1]//Expand//GetDigit; ]; (************************************************************************) (* Define GetBC[k] *) (* Get B.C. automatically based on the definition of BC[n,z,u,lambda] *) (************************************************************************) GetBC[i_,k_] := Module[{temp,j,phi,q,LAMB}, phi = Sum[u[j]*q^j,{j,0,k}]; LAMB = Sum[lambda[j]*q^j,{j,0,k}]; temp[0] = BC[i,z,phi,LAMB]; temp[1] = Series[temp[0],{q,0,k}]//Normal; temp[2] = Coefficient[temp[1],q^k]; temp[2]//Expand//GetDigit ]; (************************************************************************) (* Define functions chi[m] *) (************************************************************************) chi[m_] := If[m <= 1, 0, 1]; (************************************************************************) (* Define truncated[z^n_.] *) (************************************************************************) truncated[a_] := a /; FreeQ[a,z]; truncated[a_*f_] := a*truncated[f] /; FreeQ[a,z]; truncated[p_Plus] := Map[truncated,p]; truncated[z^n_.] := If[n > NtermMax, 0, z^n]; (************************************************************************) (* Define GetReal[] and GetDigit[] *) (************************************************************************) Naccu = 100; GetReal[c_] := N[IntegerPart[c*10^Naccu]/10^Naccu, Naccu] /; NumberQ[c]; GetDigit[c_] := N[IntegerPart[c*10^Naccu]/10^Naccu, Naccu] /; NumberQ[c]; Default[GetDigit, 1] = 0; GetDigit[p_Plus] := Map[GetDigit, p]; GetDigit[c_.*f_] := GetReal[c]*f /; NumberQ[c]; GetDigit[c_.] := 0 /; NumberQ[c] && Abs[c] < 10^(-Naccu+1); (************************************************************************) (* Define int[f,x,x0,x1,Nintegral] *) (* Integration of f in the interval [x0,x1] by Nintegral points *) (************************************************************************) int[f_, x_, x0_, x1_, Nintegral_] := Module[{temp,dx,s,t,i,M}, M = Nintegral; dx = N[x1-x0,100]/M; temp[0] = Series[f,{x,0,2}]//Normal; temp[0] = temp[0] /. x^_. ->0 ; s = temp[0]; For[i=1,i<=M,i=i+1, t = x0 + i*dx; temp[i] = f /.x->t; s = s + temp[i]//Expand; ]; temp[M+1] = (temp[0]+temp[M])/2; (s - temp[M+1])*dx//Expand ]; (************************************************************************) (* Define ChebyApproxA[f,x,a,b,M] *) (* Approximate a function by Chebyshev polynomial *) (************************************************************************) ChebyApproxA[f_,x_,a_,b_,Ntruncated_] := Module[{temp,n,z,t,F,A}, temp[0] = f/. x-> a + (b-a)*(z+1)/2; F = temp[0] /. z -> Cos[t]; For[n=0, n<=Ntruncated, n++, A[n] = 2*int[F*Cos[n*t],t,0,Pi,100]/Pi ]; For[n = 0, n<=Ntruncated, n++, temp[0] = A[n]//Expand; A[n] = GetDigit[temp[0]]; ]; temp[1] = A[0]/2 + Sum[ A[n] * ChebyshevT[n,z] , {n,1,Ntruncated}]; temp[1] /. z-> -1 + 2*(x-a)/(b-a)//Expand ]; (************************************************************************) (* Define TrigApprox[f,x,Ntruncated,TypeBase] *) (* Approximate f[x] in [0,xR] by the hybrid-base method *) (* TypeBase = 0 : pure Fourier series is used *) (* TypeBase = 1 : hybrid-base method is used *) (* Ntruncated : number of truncated terms *) (************************************************************************) TrigApprox[f_,x_,xR_,Ntruncated_,TypeBase_]:= Module[{temp,a,b,c,y,i,j,t}, If[TypeBase == 0, temp[0] = f/.x->t]; If[TypeBase == 1, temp[0] = Series[f,{x,0,2}]//Normal; temp[1] = temp[0] /. x^_. -> 0; temp[2] = f /. x->xR; temp[3] = temp[1] + (temp[2]-temp[1])/xR*x//Expand; y = GetDigit[temp[3]]; temp[0] = f - y /. x->t; ]; If[TypeBase == 2, temp[0] = D[f,x]//Expand; temp[1] = Series[temp[0],{x,0,2}]//Normal; temp[1] = temp[1] /. x^_. -> 0; temp[2] = temp[0] /. x->xR; a = temp[1]; b = -(temp[1] + temp[2])/2/xR; temp[3] = (a*x+b*x^2)*Cos[Pi*x/xR]//Expand; y = GetDigit[temp[3]]; temp[0] = f - y /. x-> t ; ]; If[TypeBase == 0 || TypeBase == 2, For[i = 0, i <= Ntruncated, i++, c[i] = 2*int[temp[0]*Cos[i*t*Pi/xR],t,0,xR,Nintegral]/xR; ]; temp[4] = c[0]/2 + Sum[c[j]*Cos[j*x*Pi/xR],{j,1,Ntruncated}]; ]; If[TypeBase == 1, For[i = 1, i <= Ntruncated, i++, c[i] = 2*int[temp[0]*Sin[i*t*Pi/xR],t,0,xR,Nintegral]/xR; ]; temp[4] = Sum[c[j]*Sin[j*x*Pi/xR],{j,1,Ntruncated}]; ]; If[TypeBase == 0, temp[5] = GetDigit[temp[4]], temp[5] = GetDigit[temp[4] + y] ]; temp[5] ]; (************************************************************************) (* Define GetMin1D[f,x,a,b,Npoint] *) (************************************************************************) GetMin1D[f_,x_,x0_,x1_,Npoint_] := Module[ {temp,fmin,Xmin,xmin0,xmin1,dx,Num,i,j,X,s}, Num = 0; dx = (x1-x0)/Npoint //GetDigit; For[i = 0, i <= Npoint, i++, X[i] = x0 + i*dx //GetDigit; temp[i] = f /. x -> X[i]; ]; For[i = 1, i <= Npoint-1, i++, If[temp[i] < temp[i-1] && temp[i] < temp[i+1], Num = Num + 1; fmin[Num] = temp[i]; Xmin[Num] = X[i]; ]; ]; For[i = 1, i <= Num, i++, j = 0; xmin0 = Xmin[i]; Label[100]; j = j + 1; temp[0] = xmin0 - dx/5^(j-1) //GetDigit; temp[1] = xmin0 + dx/5^(j-1) //GetDigit; s = GetMin1D0[f,x,temp[0],temp[1],10]; xmin1 = s[[2]]; If[Abs[xmin1 - xmin0] < 10^(-20), Xmin[i] = xmin1; fmin[i] = f /. x-> xmin1; Goto[200]; ]; xmin0 = xmin1; Goto[100]; Label[200]; Print[" Minimum = ",fmin[i]//N, " at ", x, " = ",N[Xmin[i],20]]; ]; ]; GetMin1D0[f_,x_,x0_,x1_,Npoint_] := Module[{temp, fmin, Xmin, dx, i,z,F}, dx = (x1-x0)/Npoint // GetDigit; fmin = 10^100; Xmin = 0; For[i = 1, i <= Npoint, i++, z = x0 + i*dx // GetDigit; F = f /. x->z // GetDigit; If[F < fmin, fmin = F; Xmin = z; ]; ]; {fmin,Xmin} ]; GetMin2D[f_,x_,x0_,x1_,y_,y0_,y1_,Npoint_] := Module[{temp,fmin,Xmin,Ymin,xmin0,xmin1,ymin0,ymin1,dx,dy,Num,i,j,X,Y,s}, Num = 0; dx = (x1-x0)/Npoint //GetDigit; dy = (y1-y0)/Npoint //GetDigit; For[i = 0, i <= Npoint, i++, X[i] = x0 + i*dx //GetDigit; For[j = 0, j <= Npoint, j++, Y[j] = y0 + j*dy //GetDigit; temp[i,j] = f /. {x -> X[i], y->Y[j]}; ]; ]; For[i = 1, i <= Npoint-1, i++, For[j = 1, j <= Npoint-1, j++, If[temp[i,j] < temp[i-1,j] && temp[i,j] < temp[i+1,j] && temp[i,j] < temp[i,j-1] && temp[i,j] < temp[i,j+1], Num = Num + 1; fmin[Num] = temp[i,j]; Xmin[Num] = X[i]; Ymin[Num] = Y[j]; ]; ]; ]; For[i = 1, i <= Num, i++, j = 0; xmin0 = Xmin[i]; ymin0 = Ymin[i]; Label[100]; j = j + 1; X[0] = xmin0 - dx/5^(j-1) //GetDigit; X[1] = xmin0 + dx/5^(j-1) //GetDigit; Y[0] = ymin0 - dy/5^(j-1) //GetDigit; Y[1] = ymin0 + dy/5^(j-1) //GetDigit; s = GetMin2D0[f,x,X[0],X[1],y,Y[0],Y[1],10]; xmin1 = s[[2]]; ymin1 = s[[3]]; If[Abs[xmin1 - xmin0] < 10^(-20) && Abs[ymin1 - ymin0] < 10^(-20), Xmin[i] = xmin1; Ymin[i] = ymin1; fmin[i] = f /. {x-> xmin1, y->ymin1}; Goto[200]; ]; xmin0 = xmin1; ymin0 = ymin1; Goto[100]; Label[200]; Print[" Minimum = ",fmin[i]//N]; Print[" at ", x, " = ",N[Xmin[i],20]]; Print[" ", y, " = ",N[Ymin[i],20]]; ]; ]; GetMin2D0[f_,x_,x0_,x1_,y_,y0_,y1_,Npoint_] := Module[{temp,fmin,Xmin,Ymin,dx,dy,i,X,Y,F}, dx = (x1-x0)/Npoint // GetDigit; dy = (y1-y0)/Npoint // GetDigit; fmin = 10^10; Xmin = 10^10; Ymin = 10^10; For[i = 1, i <= Npoint, i++, X = x0 + i*dx // GetDigit; For[j = 1, j <= Npoint, j++, Y = y0 + j*dy // GetDigit; F = f /. {x->X, y->Y} // GetDigit; If[F < fmin, fmin = F; Xmin = X; Ymin = Y ]; ]; ]; {fmin,Xmin,Ymin} ]; (************************************************************************) (* Define GetErr[m] *) (* This module gives averaged squared residual *) (************************************************************************) GetErr[k_] := Module[{temp}, error[k] = f[z,U[k],Lambda[k-1]]; If[ComplexQ == 0, temp[0] = error[k]^2 ] ; If[ComplexQ == 1, temp[1] = Re[error[k]]; temp[2] = Im[error[k]]; temp[0] = temp[1]^2 + temp[2]^2; ]; If[!NumberQ[zRintegral], If[FreeQ[zR,infinity] && zR < 100, zRintegral = zR, Print["Squared residual is integrated in the interval [",zL,",b]"]; Print[" The value of b = ? "]; zRintegral = Input[]; ]; ]; temp[1] = Abs[zRintegral-zL]; Err[k] = int[temp[0],z,zL,zRintegral,Nintegral]/temp[1]; ]; (************************************************************************) (* Define hp[f_,m_,n_] *) (************************************************************************) hp[f_,m_,n_]:=Block[{k,i,df,res,q}, df[0] = f[0]; For[k = 1, k <= m+n, k++, df[k] = f[k] - f[k-1]//Expand ]; res = df[0] + Sum[df[i]*q^i,{i,1,m+n}]; Pade[res,{q,0,m,n}]/.q->1 ]; (************************************************************************) (* Main Code *) (* HAM approach without iteration *) (************************************************************************) BVPh[begin_,end_]:=Block[ {uSpecial,CC,temp,ss,w,i,j,EQ,Unknown,phi,LAM,sLeng,sMin,sMinI,RHS,base,a}, time[0] = SessionTime[]; For[k=begin,k<=end,k=k+1, If[PRN == 1, Print["k = ",k]]; If[k == 1 && NumberQ[c0], c0 = GetDigit[c0]]; If[k == 1, Clear[lambda]]; (* Get special solution *) Getdelta[k-1]; RHS = H[z]*delta[k-1]//Expand; temp[1] = Linv[RHS]; uSpecial = c0*temp[1] + chi[k]*u[k-1]//Expand; (* Get lambda[k-1] and u[k] *) temp[1] = DSolve[L[w[z]]==0,w[z],z]; temp[2] = temp[1][[1,1,2]]; Unknown = {}; For[i = 1, i <= OrderEQ, i++, base[i] = Coefficient[temp[2],C[i]]; If[FreeQ[base[i],Exp[n_.*z]], Unknown = Union[Unknown,{CC[i,k]}], If[( Abs[base[i]] /. z -> 1. ) < 1, Unknown = Union[Unknown,{CC[i,k]}], CC[i,k] = 0; base[i] = 0; ]; ]; ]; ss = Sum[base[i]*CC[i,k],{i,1,OrderEQ}]; u[k] = uSpecial + ss // Expand; For[i = 1, i <= OrderEQ, i++, eq[i] = GetBC[i,k]//Expand//GetDigit]; EQ = {}; For[i=1, i <= OrderEQ, i++, If[FreeQ[BC[i,z,w[z],a],infinity], EQ = Union[EQ,{eq[i] == 0}]; ]; ]; If[TypeEQ != 1, eq[0] = GetBC[0,k]//Expand//GetDigit; EQ = Union[EQ,{eq[0] == 0}]; Unknown = Union[Unknown,{lambda[k-1]}]; ]; ss = Solve[EQ,Unknown]; sLeng = Length[ss]; If[sLeng == 1, sNum = 1]; If[k == 1 && !NumberQ[sNum], For[j = 1, j <= sLeng, j++, temp[1] = lambda[0] /. ss[[j]]; Print[j," th solution of lambda[0] = ", temp[1]//N]; ]; Print["which solution is chosen ? "]; sNum = Max[1,Input[]]; ]; If[k == 1 && TypeEQ != 1 && NumberQ[sNum] && NumberQ[LAMBDA[nIter-1]], sMin = 10^10; For[j = 1, j <= sLeng, j++, temp[1] = lambda[0] /. ss[[j]]; temp[2] = Abs[temp[1] - LAMBDA[nIter-1]]; If[ temp[2] < sMin, sNum = j; sMin = temp[2]; ]; ]; ]; For[i = 1, i <= OrderEQ, i++, temp[1] = CC[i,k] /. ss[[sNum]]//Expand; CC[i,k] = temp[1]//GetDigit; ]; u[k] = u[k]//Expand//GetDigit; U[k] = Expand[U[k-1] + u[k]]; Uz[k] = D[U[k],z]; If[TypeEQ != 1, temp[1] = lambda[k-1] /. ss[[sNum]]//Expand; lambda[k-1] = temp[1]//GetDigit; Lambda[k-1] = Sum[lambda[i],{i,0,k-1}]//Expand//GetDigit; If[PRN == 1, Print[k-1,"th-order approx. of eigenvalue = ",Lambda[k-1]//N ]; ]; ]; (* Print results *) If[PRN == 1 && NumberQ[c0] && TypeEQ == 1, output[z,U,k] ]; If[IntegerQ[k/NgetErr] && PRN == 1, GetErr[k]; If[NumberQ[Err[k]], Print["Squared residual = ",Err[k]//N]; ]; If[Err[k] < ErrReq, Print["Congratulation: the required squared residual is satisfied ! "]; Goto[end] ]; ]; If[PRN == 1, time[k] = SessionTime[]; temp[0] = time[k]-time[0]; Print["Used CPU times = ",temp[0], " (seconds) "]; ]; ]; If[PRN == 1, Print[" Successful !"]]; Label[end]; ]; (************************************************************************) (* Main Code *) (* HAM approach with iteration *) (* begin : number of starting iteration *) (* end : number of stopping iteration *) (* OrderIter : order of iteration formula *) (************************************************************************) iter[begin_, end_, OrderIter_]:=Block[{temp,time}, PRN = 0; If[TypeL != 1, ApproxQ = 1]; For[ nIter = Max[1,begin], nIter <= end, nIter = nIter + 1, If[nIter == 1, sNum = .; If[NumberQ[c0], c0 = GetDigit[c0]]; Time[0] = SessionTime[]; V[0] = u[0]//Expand//GetDigit; Vz[0] = D[V[0],z]; ]; If[nIter > 1, u[0] = V[nIter-1]; U[0] = u[0]; ]; Print[ nIter, " th iteration:"]; If[NumberQ[c0], BVPh[1,OrderIter]]; If[!NumberQ[c0], BVPh[1,Max[3,OrderIter]]; GetErr[1]; GetErr[2]; GetErr[3]; Print[" Squared Residual of Governing Equation : "]; Print[" Red line : 1st-order approx."]; Print[" Green line : 2nd-order approx."]; Print[" Blue line : 3rd-order approx."]; LogPlot[{Abs[Err[1]],Abs[Err[2]],Abs[Err[3]]},{c0,c0L,0}, PlotStyle->{RGBColor[1,0,0],RGBColor[0,1,0],RGBColor[0,0,1]}]; Print[" Adjust the interval of c0 by choosing a better end-point "]; Print[" New end-point on left (c0 < 0) or right (c0 > 0) = ? (Input 0 to skip it) "]; c0L = Input[]; If[ c0L != 0, LogPlot[{Abs[Err[1]],Abs[Err[2]],Abs[Err[3]]},{c0,c0L,0}, PlotStyle->{RGBColor[1,0,0],RGBColor[0,1,0],RGBColor[0,0,1]}]; ]; Print[" Choose the value of convergence-control parameter c0 "]; Print[" c0 = ? ( Input 0 to stop the code) "]; temp[0] = Input[]; c0 = GetDigit[temp[0]]; Print[" c0 = ",c0//N]; If[ c0 == 0, c0 =.; Print["Stop!"]; Goto[stop]]; ]; temp[0] = U[OrderIter]//Expand; If[PolynomialQ[temp[0],z], temp[0] = truncated[temp[0]]]; V[nIter] = GetDigit[temp[0]]; Vz[nIter] = D[V[nIter],z]; If[TypeEQ == 1, output[z,V,nIter], LAMBDA[nIter] = Lambda[OrderIter-1]//Expand; Print[" Eigenvalue = ",LAMBDA[nIter]//N ]; ]; GetErr[OrderIter]; ERR[nIter] = Err[OrderIter]; Print[" Squared Residual = ",ERR[nIter]//N]; If[ERR[nIter] < ErrReq, Print[" STOP : Required accuracy is satisfied ! "]; Goto[stop]; ]; If[ nIter > 1 && ERR[nIter] > ERR[nIter-1] && Err[nIter] < 10^(-7), Print[" Squared residual does NOT decrease !"]; Goto[stop]; ]; Time[nIter] = SessionTime[]; Temp[-1] = Time[nIter]-Time[0]; Print[" Used CPU times = ", Temp[-1], " (seconds) "]; If[IntegerQ[nIter/Nupdate] && Abs[c0]< 1/2, c0 =.]; ]; Print[" End of iteration ! "]; Label[stop]; PRN = 1; ]; (************************************************************************) (* Print equation, boundary conditions and control parameters *) (************************************************************************) PrintInput[s_] := Module[{}, Print["-----------------------------------------------------------------"]; Print["The values of control parameters: "]; Print[" Nupdate = ", Nupdate]; Print[" Ntruncated = ", Ntruncated]; Print[" NtermMax = ", NtermMax]; Print[" Nintegral = ", Nintegral]; Print[" ErrReq = ", ErrReq//N]; Print[" NgetErr = ", NgetErr]; Print[" ComplexQ = ", ComplexQ]; Print[" PRN = ", PRN]; Print[" c0L = ", c0L]; Print[" zL = ", zL]; Print["-----------------------------------------------------------------"]; Print[" Governing Equation : ", f[z,s,lambda], " = 0 "]; Print[" Interval of z : ",zL, " < z < ", zR]; If[TypeEQ != 1, Print[" ",0,"th Boundary Condition : ", BC[0, z, s, lambda]," = 0"]; ]; For[i=1, i<= OrderEQ, i++, Print[" ",i,"th Boundary Condition : ", BC[i, z, s, lambda]," = 0"]; ]; Print["-----------------------------------------------------------------"]; If[ApproxQ == 0, Print[" Delta[k] is NOT approximated by other base functions ! "], If[TypeL == 1, Print[" Base function : Chebyshev "], If[TypeBase == 1, Print[" Base function : hybrid-base with sine and polynomial"], Print[" Base function : hybrid-base with cosine and polynomial "] ]; ]; ]; Print[" Auxiliary linear operator L[u] = ", L[s] ]; If[TypeL != 1, Print[" kappa = ", kappa]]; Print[" Initial guess u[0] = ", u[0] ]; Print[" Auxiliary function H[z] = ", H[z]]; If[NumberQ[c0],Print[" Convergence-control parameter c0 = ", c0]]; Print["-----------------------------------------------------------------"]; ];