(************************************************************************)
(*             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  *) 
(************************************************************************)
<<Calculus`Pade`;
<<Graphics`Graphics`;


(************************************************************************)
(*                      Define initial guess of w(x)                    *)
(************************************************************************)
w[0] = 0;
U[0] = A + w[0];

(************************************************************************)
(*                     Define the function chi[k]                       *)
(************************************************************************)
chi[k_]:=If[ k <= 1, 0, 1];

(************************************************************************)
(*              Define the the auxiliary linear operator L              *)
(************************************************************************)
L[f_] := Module[{temp},
Expand[c2*D[f,{x,1},{y,1}]/x/y + c4*D[f,{x,2},{y,2}]] 
];


(************************************************************************)
(*    Define the inverse operator of the auxiliary linear operator L    *)
(************************************************************************)
Linv[x^m_. * y^n_.] := x^(m+2)*y^(n+2)/(m+2)/(n+2)/(c2+c4*(n+1)*(m+1));
Linv[x^m_.] := x^(m+2) * y^2/(m+2)/2/(c2+c4*(m+1));
Linv[y^n_.] := x^2 * y^(n+2)/(n+2)/2/(c2+c4*(n+1));
Linv[c_]    := c*x^2*y^2/4/(c2+c4) /;   FreeQ[c,x] && FreeQ[c,y];


(************************************************************************)
(*              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,x] && FreeQ[c,y];

(************************************************************************)
(*                        Define Getdelta[k]                            *)
(************************************************************************)
Getdelta[k_]:=Module[{temp,n},
temp[1] = D[w[k],{x,2}] + D[w[k],{y,2}];
temp[2] = Sum[lambda[n]*Dexp[k-n],{n,0,k}];
delta[k]    = Expand[ temp[1] + temp[2] ];
];

(************************************************************************)
(*                             Define GetDexp[n]                       *)
(************************************************************************)
GetDexp[n_]:=Module[{j},
If[n == 0, Dexp[0] = Exp[w[0]] ];
If[n  > 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]]