(************************************************************************) (* Mathematica code for Example 2.2 *) (* Solve a nonlinear oscillation governed by: *) (* x'' + lambda* x + epsilon*x^3 = 0 *) (* subject to the initial conditions: *) (* x(0) = x*, x'(0) = 0 *) (* VERSION A: without iteration *) (* *) (* Dr. Shijun Liao *) (* Shanghai Jiao Tong University *) (* China *) (* June, 2010 *) (************************************************************************) < 0 && lambda >= 0, beta = 0]; If[epsilon > 0 && lambda < 0 && xstart > 0, beta = N[Sqrt[Abs[lambda/epsilon]],100]]; If[epsilon > 0 && lambda < 0 && xstart < 0, beta =-N[Sqrt[Abs[lambda/epsilon]],100]]; If[epsilon < 0 && lambda >= 0 && Abs[xstart] < Sqrt[Abs[lambda/epsilon]], beta = 0]; If[!NumberQ[beta//N],Print[" There is no periodic solution for these input data !"]]; If[xOptimal == 1, beta = . ]; (************************************************************************) (* Define initial guess of x(tau) *) (************************************************************************) x[0] = beta + (xstart - beta)*Cos[t]; X[0] = x[0]; u[0] = X[0] /. t->Sqrt[GAMMA[0]]*t ; (************************************************************************) (* Define the function chi[k] *) (************************************************************************) chi[k_] := If[ k <= 1, 0, 1 ]; (************************************************************************) (* Define the inverse operator of the auxiliary linear operator L *) (************************************************************************) Linv[Cos[m_*t]] := (-1)^(KAPPA+1)*Cos[m*t]/(1-m^(2*KAPPA)); (************************************************************************) (* The linear 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,t]; Linv[c_] := (-1)^(KAPPA+1)*c /; FreeQ[c,t]; (************************************************************************) (* Define GetDelta[k] *) (* Calculate the right-hand side term delta[k] according to (2.55) *) (************************************************************************) GetDelta[k_]:=Module[{temp}, temp[1] = xttgamma[k]//TrigReduce; temp[2] = xxx[k]//TrigReduce; temp[3] = temp[1] + lambda*x[k] + epsilon * temp[2]; Delta[k]= Expand[temp[3]]; ]; (************************************************************************) (* Define GetxAll *) (************************************************************************) GetxAll[k_]:=Module[{temp}, xtt[k] = Expand[D[x[k],{t,2}]]; temp = Sum[x[j]*x[k-j],{j,0,k}]//Expand; xx[k] = TrigReduce[temp]; temp = Sum[x[j]*xx[k-j],{j,0,k}]//Expand; xxx[k] = TrigReduce[temp]; xttgamma[k] = Sum[gamma[j]*xtt[k-j],{j,0,k}]//Expand; ]; (************************************************************************) (* Define Getlambda *) (* This module gets lambda[k-1] *) (************************************************************************) Getgamma[k_]:=Module[{temp,eq}, temp[1] = Delta[k]//TrigReduce; temp[2] = Expand[temp[1]]; eq = Coefficient[temp[2],Cos[t]]; temp[3] = Solve[eq == 0, gamma[k]]; gamma[k] = temp[3][[1,1,2]]; ]; (************************************************************************) (* Define Getx *) (* Get a solution of high-order deformation equation *) (************************************************************************) Getx[k_]:=Module[{temp}, temp[1] = xSpecial + chi[k]*x[k-1]; temp[2] = temp[1] /. t->0; x[k] = temp[1] - temp[2]*Cos[t]; ]; (************************************************************************) (* Define GetErr *) (* Get residual error square of governing equation *) (************************************************************************) GetErr[k_]:=Module[{temp,Xtt,error,delt,tt,Npoint,sum,i}, Xtt = D[X[k],{t,2}]; error = GAMMA[k]*Xtt + lambda * X[k] + epsilon * X[k]^3; Npoint = 50; sum = 0; delt = N[2*Pi/Npoint,24]; For[ i = 0, i <= Npoint, i=i+1, tt = i*delt; temp = error/.t -> tt //Expand; sum = sum + temp^2 //Expand; ]; Err[k] = sum/(Npoint+1) ]; (************************************************************************) (* Main Code *) (************************************************************************) ham[m0_,m1_]:=Module[{temp,k,j,zz}, For[k=Max[1,m0], k<=m1, k=k+1, Print[" k = ",k]; GetxAll[k-1]; GetDelta[k-1]; RHS[k] = c0*Delta[k-1]//Expand; Getgamma[k-1]; GAMMA[k-1] = Expand[Sum[gamma[j],{j,0,k-1}]]; T[k-1] = 2*Pi/Sqrt[GAMMA[k-1]]; If[NumberQ[GAMMA[k-1]], Print[k-1,"th approx. of gamma = ",N[GAMMA[k-1],10], " variation = ",gamma[k-1]//N ]; Print[k-1,"th approx. of T = ",N[T[k-1],10] ]; ]; If[k == 1 && xOptimal == 1, GetErr[0]; If[ xstart > 0, temp = Minimize[{Err[0], beta > xstart},beta] ]; If[ xstart < 0, temp = Minimize[{Err[0], beta < xstart},beta] ]; zz = beta/.temp[[2]]; beta = N[IntegerPart[zz*10000]/10000,100]; Print[" Optimal initial approx. is used with beta = ", beta//N]; ]; temp = RHS[k]/.Cos[t]->0//Expand; xSpecial = Linv[temp]; Getx[k]; X[k] = X[k-1] + x[k]//Expand; u[k] = X[k] /. t-> Sqrt[GAMMA[k-1]]*t; ]; Print["Successful !"]; ]; (************************************************************************) (* Print the physical and auxiliary parameters *) (************************************************************************) Print[" lambda = ",lambda]; Print[" epsilon = ",epsilon]; Print[" x* = ",xstart]; Print[" beta = ",beta]; Print[" c0 = ",c0]; Print[" KAPPA = ",KAPPA]; Print[" xOptimal = ",xOptimal]; (* Gain 10th-order hootopy-approximation *) ham[1,20];