(************************************************************************) (* Mathematica code for 2D Gelfand equation *) (* *) (* For given A, it gives such an eigenvalue lambda and a normalized *) (* eigenfunction w(x,y) satisfying: *) (* w_{xx} + w_{yy} + lambda * Exp[w(x)] = 0 *) (* subject to the boundary conditions: *) (* w(1,y) = w(-1,y) = w(x,1) = w(x,-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}; temp = Solve[eq == 0, lambda[k-1]]; lambda[k-1] = temp[[1,1,2]]//Expand; ]; (************************************************************************) (* Define GetwSpecial *) (* This module gets a special solution of Eq. (8.23) *) (************************************************************************) GetwSpecial[k_]:=Module[{temp}, temp[0] = Expand[RHS[k]]; temp[1] = Linv[temp[0]]; temp[2] = temp[1] + chi[k]*w[k-1]//Simplify; wSpecial = temp[2]//Expand; ]; (************************************************************************) (* Define Getw *) (* This module gets a solution of Eq. (8.23) and (8.24) *) (************************************************************************) Getw[k_]:=Module[{temp,alpha}, alpha = (c0+chi[k])*w[k-1]+c0*(1-chi[k])*A; temp[1] = wSpecial /. y->1; temp[2] = wSpecial /. x->1; temp[3] = wSpecial /. {x->1, y->1}; temp[4] = alpha /. x->1; temp[5] = alpha /. y->1; temp[6] = alpha /. {x->1,y->1}; temp[7] = wSpecial-temp[1]-temp[2]+temp[3]+temp[4]+temp[5]-temp[6]; w[k] = Simplify[temp[7]]//Expand; ]; (************************************************************************) (* Define GetErr[k] *) (* This module gets a solution of Eq. (8.23) and (8.24) *) (************************************************************************) GetErr[k_]:=Module[{temp,sum,dx,dy,Num,i,j,X,Y}, err[k] = D[U[k],{x,2}] + D[U[k],{y,2}] + LAMBDA[k]*Exp[U[k]]; Nx = 10; Ny = 10; dx = N[1/Nx,100]; dy = N[1/Ny,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; temp = err[k]^2 /. {x->X, y->Y}; sum = sum + temp; Num = Num + 1; ]; ]; Err[k] = sum/Num; If[NumberQ[Err[k]], Print["Squared Residual of G.E. = ", Err[k]//N]] ]; (************************************************************************) (* define Body1[A] *) (* Boyd's one-point analytic approximation *) (************************************************************************) Boyd1[A_] := 3.2*A*Exp[-0.64*A]; (************************************************************************) (* define Body3[A] *) (* Boyd's three-point analytic approximation *) (************************************************************************) Boyd3[A_] := Module[{temp,B,C,G}, G = 0.2763 + Exp[0.463*A] + 0.0483*Exp[-0.209*A]; B = A*( 0.829-0.566*Exp[0.463*A]-0.0787*Exp[-0.209*A])/G; C = A*(-1.934+0.514*Exp[0.463*A]+1.9750*Exp[-0.209*A])/G; temp[1] = 2.667*A + 4.830*B + 0.127*C; temp[2] = 0.381*A + 0.254*B + 0.018*C; temp[1]*Exp[-temp[2]] ]; (************************************************************************) (* 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]//Simplify; Getlambda[k]; Lambda[k-1] = Sum[lambda[j],{j,0,k-1}]//Expand; LAMBDA[k-1] = Lambda[k-1]*Exp[-A]; Print[k-1,"th approximation of LAMBDA = ",N[LAMBDA[k-1],10] ]; ]; Print["Successful !"]; ]; (************************************************************************) (* Define the parameters *) (************************************************************************) c0 = -1; A = 1; c2 = 0; c4 = -1/4; (************************************************************************) (* Print the parameters, initial guessesand auxiliary function *) (************************************************************************) Print[" A = ",A]; Print[" c0 = ",c0]; Print[" c2 = ",c2]; Print[" c4 = ",c4]; (* Get the 20th-order approximation of the eigenvalue LAMBDA *) ham[1,21]; (* Get squared residual of up to 20th-order approx. *) For[k=5, k<=20, k=k+5, Print["k = ",k]; GetErr[k]]