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

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

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


(************************************************************************)
(*    Define the inverse operator of the auxiliary linear operator L    *)
(************************************************************************)
Linv[x^m_ * y^n_ * z^k_] := x^(m+2)*y^(n+2)*z^(k+2)/(m+2)/(n+2)/(k+2)/(c3+c6*(n+1)*(m+1)*(k+1));
Linv[y^n_ * z^k_] := x^2*y^(n+2)*z^(k+2)/2/(n+2)/(k+2)/(c3+c6*(n+1)*(k+1));
Linv[x^m_ * z^k_] := x^(m+2)*y^2*z^(k+2)/(m+2)/2/(k+2)/(c3+c6*(m+1)*(k+1));
Linv[x^m_ * y^n_] := x^(m+2)*y^(n+2)*z^2/(m+2)/(n+2)/2/(c3+c6*(n+1)*(m+1));
Linv[x^m_] := x^(m+2)*y^2*z^2/(m+2)/2/2/(c3+c6*(m+1));
Linv[y^n_] := x^2*y^(n+2)*z^2/2/(n+2)/2/(c3+c6*(n+1));
Linv[z^k_] := x^2*y^2*z^(k+2)/2/2/(k+2)/(c3+c6*(k+1));
Linv[c_]   := c*x^2*y^2*z^2/8/(c3+c6) /;   FreeQ[c,x] && FreeQ[c,y] && FreeQ[c,z];


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

(************************************************************************)
(*                            Define GetR[k]                            *)
(*                  Get the term R[k] defined by (8.25)                 *)
(************************************************************************)
Getdelta[k_]:=Module[{temp,n},
temp[1]  = D[w[k],{x,2}] + D[w[k],{y,2}] + D[w[k],{z,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[{},
If[n == 0, Dexp[0] = 1 ];
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, 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]]