(************************************************************************) (* Mathematica code for 3D Gelfand equation *) (* *) (* For given A, it gives such an eigenvalue lambda and a normalized *) (* eigenfunction w(x,y) satisfying: *) (* w_{xx} + w_{yy} + w_{zz} + lambda * Exp[w(x)] = 0 *) (* subject to the boundary conditions: *) (* w(1,y,z) = w(-1,y,z) = w(x,1,z) = w(x,-1,z) = -A, *) (* w(x,y,1) = w(x,y,-1) = -A, w(0,0) = 0. *) (* *) (* Dr. Shijun Liao *) (* Shanghai Jiao Tong University *) (* China *) (* August, 2010 *) (************************************************************************) < 0, Dexp[n] = Sum[(1-j/n)*Dexp[j]*w[n-j],{j,0,n-1}]//Expand ]; ]; (************************************************************************) (* Define Getlambda *) (* This module gets lambda[k-1] *) (************************************************************************) Getlambda[k_]:=Module[{eq,temp}, eq = w[k]-(c0+chi[k])*w[k-1] /. {x->0, y->0, z->0}; temp = Solve[eq == 0, lambda[k-1]]; lambda[k-1] = temp[[1,1,2]]//Expand; ]; (************************************************************************) (* Define GetwSpecial *) (* This module gets a special solution *) (************************************************************************) GetwSpecial[k_]:=Module[{temp}, temp[0] = Expand[RHS[k]]; temp[1] = Linv[temp[0]]; temp[2] = temp[1] + chi[k]*w[k-1]; wSpecial = temp[2]//Expand; ]; (************************************************************************) (* Define Getw *) (* This module gets a solution of high-order deformation eq *) (************************************************************************) Getw[k_]:=Module[{temp,mu}, mu = (c0+chi[k])*w[k-1] + c0*(1-chi[k])*A; temp[1] = wSpecial /. x->1; temp[2] = wSpecial /. y->1; temp[3] = wSpecial /. z->1; temp[4] = wSpecial /. {x->1, y->1}; temp[5] = wSpecial /. {x->1, z->1}; temp[6] = wSpecial /. {y->1, z->1}; temp[7] = wSpecial /. {x->1, y->1, z->1}; temp[8] = wSpecial-temp[1]-temp[2]-temp[3]+temp[4]+temp[5]+temp[6]-temp[7]; temp[1] = mu /. x->1; temp[2] = mu /. y->1; temp[3] = mu /. z->1; temp[4] = mu /. {x->1, y->1}; temp[5] = mu /. {x->1, z->1}; temp[6] = mu /. {y->1, z->1}; temp[7] = mu /. {x->1, y->1, z->1}; temp[9] = temp[8]+temp[1]+temp[2]+temp[3]-temp[4]-temp[5]-temp[6]+temp[7]; w[k] = Expand[temp[9]]; ]; (************************************************************************) (* Define GetErr[k] *) (************************************************************************) GetErr[k_]:=Module[{temp,sum,dx,dy,dz,Num,i,j,m,X,Y,Z,Nx,Ny,Nz}, err[k] = D[U[k],{x,2}] + D[U[k],{y,2}] + D[U[k],{z,2}] + LAMBDA[k]*Exp[U[k]]; Nx = 10; Ny = 10; Nz = 10; dx = N[1/Nx,100]; dy = N[1/Ny,100]; dz = N[1/Nz,100]; sum = 0; Num = 0; For[i = 0, i <= Nx-1, i++, X = i*dx; For[j = 0, j <= Ny - 1, j++, Y = j*dy; For[m = 1, m <= Nz-1, m++, Z = m*dz; temp = err[k]^2 /. {x->X, y->Y,z->Z}; sum = sum + temp; Num = Num + 1; ]; ]; ]; Err[k] = sum/Num; If[NumberQ[Err[k]], Print["Squared Residual of G.E. = ", Err[k]//N]] ]; (************************************************************************) (* Define HP[F,m,n] *) (* This module gives [m,n] homotopy-Pade approximation of the series *) (* F[k] = sum[f[i],{i,0,k}] *) (* For details, please see section 2.3.7 on page 38 to 41 *) (************************************************************************) hp[F_,m_,n_]:=Block[{i,k,dF,temp,q}, dF[0] = F[0]; For[k = 1, k <= m+n, k = k+1, dF[k] = Expand[F[k] - F[k-1]] ]; temp = dF[0]+Sum[dF[i]*q^i,{i,1,m+n}]; Pade[temp,{q,0,m,n}]/.q->1 ]; (************************************************************************) (* Main Code *) (************************************************************************) ham[m0_,m1_]:=Module[{temp,k,j}, For[k=Max[1,m0], k<=m1, k=k+1, Print[" k = ",k]; GetDexp[k-1]; Getdelta[k-1]; RHS[k] = c0 * delta[k-1]//Expand; GetwSpecial[k]; Getw[k]; U[k] = U[k-1] + w[k]//Expand; Getlambda[k]; Lambda[k-1] = Expand[Sum[lambda[j],{j,0,k-1}]]; LAMBDA[k-1] = Lambda[k-1]*Exp[-A]; Print[k-1,"th approximation of lambda = ",N[LAMBDA[k-1],10] ]; ]; Print["Successful !"]; ]; (************************************************************************) (* Physical and control input data *) (************************************************************************) A = 1; c0 = -1; c3 = 0; c6 = 1/8; (************************************************************************) (* Print input data *) (************************************************************************) Print[" A = ",A]; Print[" c0 = ",c0]; Print[" c3 = ",c3]; Print[" c6 = ",c6]; (* Get the 10th-order approximation of the eigenvalue LAMBDA *) ham[1,11]; (* Get squared residual of up to 10th-order approx. *) For[k=2, k<=10, k=k+2, Print["k = ",k]; GetErr[k]]