(************************************************************************) (* 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 B: with 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[" No periodic solution for the input data !"]]; 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]] := Cos[m*t]/(1-m^2); (************************************************************************) (* 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_] := c /; FreeQ[c,t]; (************************************************************************) (* Define GetDelta[k] *) (************************************************************************) GetDelta[k_]:=Module[{temp,n}, temp[1] = xttgamma[k]//TrigReduce; temp[2] = xxx[k]//TrigReduce; temp[3] = temp[1] + lambda*x[k] + epsilon * temp[2]; Delta[k] = temp[3]//Expand; ]; (************************************************************************) (* 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[iter_] := Module[{temp, Xtt, error, delt, tt, Npoint, sum,i}, Xtt = D[X[iter],{t,2}]; error = GAMMA[iter]*Xtt + lambda * X[iter] + epsilon * X[iter]^3; Npoint = 50; sum = 0; delt = N[2*Pi/Npoint,100]; For[ i = 0, i <= Npoint, i=i+1, tt = i*delt; temp = error/.t -> tt //Expand; sum = sum + temp^2 //Expand; ]; Err[iter] = sum/(Npoint+1); ]; (************************************************************************) (* Define Truncation *) (* Delete the high-frequency terms of R *) (************************************************************************) Truncation[k_,Nterms_] := Module[{temp,coef,i,j}, temp[1] = RHS[k] //TrigReduce; temp[2] = temp[1] //Expand; coef[0] = temp[2] /. Cos[_] -> 0; For[i = 2, i <= Nterms, i = i + 1, coef[i] = Coefficient[temp[2],Cos[i*t]]; ]; temp[3] = coef[0] + Sum[coef[j]*Cos[j*t],{j,2,Nterms}]; RHS[k] = temp[3]//Expand; ]; (************************************************************************) (* Main Code *) (************************************************************************) ham[m0_,m1_]:=Module[{temp,k,j,zz}, For[iter = Max[1,m0], iter <= m1, iter = iter + 1, For[k = 1, k <= IterOrder, k = k+1, If[k == 1, Clear[gamma]]; GetxAll[k-1]; GetDelta[k-1]; RHS[k] = c0*Delta[k-1]//Expand; Getgamma[k-1]; Truncation[k, Nterms]; xSpecial = Linv[RHS[k]]; Getx[k]; ]; X[iter] = Sum[x[i],{i,0,IterOrder}]; GAMMA[iter] = Expand[Sum[gamma[j],{j,0,IterOrder-1}]]; variation = GAMMA[iter]-chi[iter]*GAMMA[iter-1]//N; T[iter] = 2*Pi/Sqrt[GAMMA[iter]]; u[iter] = X[IterOrder] /. t-> Sqrt[GAMMA[iter]]*t; GetErr[iter]; If[ iter == 1 , c0 = .; temp = Minimize[{Err[1],c0 < 0}, c0]; zz = c0/.temp[[2]]; c0 = N[IntegerPart[zz*10000]/10000,100]; Print["The optimal value of c0 = ",c0//N]; ]; Print[ iter, "th iteration: error = ", Err[iter]//N]; Print[ " gamma = ",N[GAMMA[iter],24], " variation = ",variation]; If[ Err[iter] < 10^(-30) || Abs[variation] < 10^(-30), Print[" Desired accuracy is arrived: iteration stops !"]; Goto[end]; ]; x[0] = X[iter]; Clear[gamma]; ]; Label[end]; Print[" OK !"]; ]; (************************************************************************) (* Print the physical and auxiliary parameters *) (************************************************************************) Print[" lambda = ",lambda]; Print[" epsilon = ",epsilon]; Print[" x* = ",xstart]; Print["----------------------------------------------------------------"]; Print[" Oder of iteration formula = ",IterOrder]; Print[" N, the number of terms in delta = ",Nterms]; (* Gain 10 homotopy-iteration approximation *) ham[1,10];