(************************************************************************) (* Mathematica code for nonlinear wave-current interaction *) (* *) (* The interaction between a highly nonlinear traveling wave and a *) (* non-uniform shear current is solved analytically by the HAM. *) (* *) (* Jun Cheng & Shijun Liao *) (* Shanghai Jiao Tong University *) (* China *) (* August, 2010 *) (************************************************************************) <0//ComplexExpand; ccc[f_] := Select[f, FreeQ[#, x] &]; (************************************************************************) (* Define Gety[k] *) (************************************************************************) Gety[k_] := Module[{temp,sol,p0,p1,plist,clist,dlist,c1,property}, temp[0] = therest[ySpecial,k]; temp[1] = temp[0]//Expand; p0 = ccc[temp[1]]; p1 = Coefficient[temp[1],Cos[x]]; temp[2] = temp[1] - p0 - p1*Cos[x]; sol = Solve[{p0 == 0,p1 == 0},{gamma[k-1],mu[k-1]}]; gamma[k-1] = gamma[k-1]/.sol[[1]]//Expand; mu[k-1] = mu[k-1]/.sol[[1]]//Expand; plist = Table[Coefficient[temp[2],Cos[i*x]],{i,2,2k+1}]//Simplify; clist = Table[plist[[i-1]]/(1-i),{i,2,2k+1}]; dlist = Table[clist[[i-1]]*(1-(-1)^i),{i,2,2k+1}]; temp[3] = Apply[Plus,dlist]; temp[4] = ySpecial /.{z->0,x->0 }; temp[5] = ySpecial /.{z->0,x->Pi}; temp[6] = Table[Exp[-i*z]*(Exp[i*I*x]+Exp[-i*I*x])/2,{i,2,2k+1}]; c1 = -1/2*(temp[4]-temp[5]+temp[3]); property = Dot[temp[6],clist] + c1*Exp[-z]*(Exp[I*x] + Exp[-I*x])/2; y[k] = Expand[ySpecial + chi[k]*y[k-1] + property]; ]; (************************************************************************) (* Define GetYs[k] *) (************************************************************************) GetYs[k_]:=Module[{temp}, temp[0] = Y[k] /. z->0; temp[1] = Integrate[temp[0], {x,-Pi,Pi}]/2/Pi; Ys[k] = temp[0] - temp[1]//Expand; ]; (************************************************************************) (* Define GetErrB[k] *) (************************************************************************) GetErrB[k_]:=Module[{temp,i,Yx,Yz,Nx,Np,dx,xx}, Yx = D[Y[k],x]; Yz = D[Y[k],z]; temp[1] = MU[k]*(1+Yx^2) + 2*(Y[k]-GAMMA[k])*Yz^2 /. z->0; temp[2] = temp[1]^2; Nx = 20; dx = N[Pi/Nx,100]; sum = 0; Np = 0; For[i = 0, i <= Nx, i++, xx = i*dx; sum = sum + temp[2] /. x->xx; Np = Np + 1; ]; ErrB[k] = sum/Np; If[NumberQ[ErrB[k]], ErrB[k] = Re[ErrB[k]]; Print["k = ",k," Squared Residual of B.C. = ",ErrB[k]//N] ]; ]; (************************************************************************) (* Define GetErr[k] *) (************************************************************************) GetErr[k_]:=Module[{temp,i,j,sum,xx,zz,Yx,Yz,Yxx,Yxz,Yzz,Nx,Nz,Np}, Yx = D[Y[k],x]; Yz = D[Y[k],z]; Yxx = D[Yx,x]; Yxz = D[Yx,z]; Yzz = D[Yz,z]; temp[1] = Yxx*Yz^2 - 2*Yx*Yz*Yxz + (1+Yx^2)*Yzz - Yz^3*Omega; temp[2] = temp[1]^2; Nx = 10; Nz = 10; dx = N[ Pi/Nx,100]; dz = N[2*Pi/Nz,100]; sum = 0; Np = 0; For[i = 0, i <= Nx, i++, For[j = 0, j <= Nz,j++, xx = i*dx; zz = j*dz; sum = sum + temp[2]/.{x->xx,z->zz}; Np = Np + 1; ]; ]; Err[k] = sum/Np; If[NumberQ[Err[k]], Err[k] = Re[Err[k]]; Print["k = ",k," Squared Residual of G.E. = ",Err[k]//N] ]; ]; (************************************************************************) (* 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[m0_,m1_]:=Module[{temp,k,n}, For[k=Max[1,m0],k<=m1,k=k+1,Print[" k = ",k]; GetAll[k-1]; GetRHS[k]; GetRHSb[k]; GetySpecial[k]; Gety[k]; Y[k] = Y[k-1] + y[k]//ComplexExpand; If[k>1, MU[k-1] = MU[k-2] + mu[k-1]//Expand; GAMMA[k-1] = GAMMA[k-2] + gamma[k-1]//Expand; ]; Print[" mu = ",MU[k-1]//N," variation = ",mu[k-1]//N]; Print[" gamma = ",GAMMA[k-1]//N," variation = ",gamma[k-1]//N]; ]; Print[" Sucessful ! "]; ]; (************************************************************************) (* Initial guess of u[0] and related functions *) (************************************************************************) Y[0] = y[0]; GAMMA[0] = gamma[0]; MU[0] = mu[0]; ERR[0] = ComplexExpand[(Y[0] - GAMMA[0]) /.x->0/.z->0]; Omega = epsilon* Exp[-z]; (* Physical and control parameters *) epsilon = 1/5; H = 3/10; c0 =-11/20; (* Print input data *) Print["epsilon = ",epsilon]; Print[" H = ",H]; Print[" c0 = ",c0]; (* Gain the 10th-order approximation *) ham[1,11]; (* Gain squared residual of governing equation *) For[k=2, k<=10, k=k+2, GetErr[k]]; (* Gain squared residual of boundary condition *) For[k=2, k<=10, k=k+2, GetErrB[k]];