Trying to use Matlab to find Numerical Solution to $u''(x)+e^{u(x)}=0, u(0)=0, u(1)=0$ - Newton's method

Question

I am trying to use Matlab to find Numerical Solution to $u''(x)+e^{u(x)}=0, u(0)=0, u(1)=0$. I can't get my code to work out. I was hoping you could help me. I have added citations to my code to explain what is going on. My confusion is that I am getting a lot of different "solution" graphs output. They change each time I run it. This doesn't seem right to me. Also, the boundary conditions are not portrayed on the graphed "solution".

For $G=\frac{1}{h^2}(u_{j-1}-2u_j+u_{j+1})+N(u_j)$ where $N(u_j)$ is the $e^{u(x)}$ nonlinear part, and for J being the Jacobian, I am trying to follow this procedure: Solve for delta, $$J(U^{(k)})\delta = -G(U^{(k)})$$ Then set $$U^{(k+1)}=U^{(k)}+\delta$$ Test for accuracy and repeat if necessary epsilon isn't reached. Thank you for any help you can give.

Note: I think my initial guess should be this due to boundary conditions. Uk = [0 10*randn(1, m2-2) 0]' But this doesn't converge.

This Uk leads to convergence at least... Uk = rand(m2,1);

ax = 0;
bx = 1;
sigma = 0;   % Dirichlet boundary condition at ax
beta = 0;     % Dirichlet boundary condtion at bx
lambda = 1;
f = @(x) 0*x;  % right hand side function
epsilon = 0.01;
count = 1;
m1vals = 40;
m2= m1vals+1;
m1 = m1vals;
m = m1 - 1;                 % number of interior grid points
%Uk = [0 10*randn(1, m2-2) 0]'
%This Uk leads to convergence at least.
Uk = rand(m2,1);
Uk1=Uk;
h=0.01;
N = @(u) lambda*exp(u);
% set grid points:  
gridchoice = 'uniform';   
x = xgrid(ax,bx,m,gridchoice); 

while norm(Uk1-Uk) > epsilon || count ==1
    norm(Uk1-Uk) %check norm and print
    if count ~= 1 %only after first iteration
        Uk=Uk1;
    end
    %Part of J with -2 on main and 1 on off digaonals
    Jpart = (1/h^2)*(diag(-2*ones(m2,1),0)+diag(ones(m2-1,1),1)+diag(ones(m2-1,1),-1));
    %Create N' part of Jacobian
    %addJ is an estimate of u'(x) since we need u'(x)*lambda*exp(u(x))
    addJ = spalloc(m2,m2,3*m2);
    addJ(1,1:3) = fdcoeffF(1, x(1), x(1:3)); 
    for i=2:m1
         addJ(i,i-1:i+1) = fdcoeffF(1, x(i), x((i-1):(i+1)));
    end
    addJ(m2,m:m2) = fdcoeffF(1,x(m2),x(m:m2)); 
    %computes u'(x)*lambda*exp(u(x)) at Uk.
    J2 = diag(addJ*Uk,0)*N(Uk);   %addJ*Uk is an estimate of u'(x) at Uk points.
    J=Jpart+diag(J2,0); %put all together

      % set up matrix G = u''+N
      Gpart = spalloc(m2,m2,3*m2);   % initialize to zero matrix
      % first row for Dirichlet BC at ax:
      Gpart(1,1:3) = fdcoeffF(0, x(1), x(1:3)); 
      % interior rows:
      for i=2:m1
         Gpart(i,i-1:i+1) = fdcoeffF(2, x(i), x((i-1):(i+1)));
      end
      % last row for Dirichlet BC at bx:
      Gpart(m2,m:m2) = fdcoeffF(0,x(m2),x(m:m2)); 
      %Set up for J(Uk)delta=-G(Uk) 
      Gall=Gpart+diag(N(Uk),0); 
      Gans=-Gall*Uk; %evaluate G(Uk)
Delta = J\Gans;
Uk1 = Uk + Delta;
count = count + 1;  
end



  clf
  plot(x,Uk1,'o')  % plot computed solution
  title(sprintf('Computed solution with %i grid points',m2));

Answer 1

I had trouble following your code, but it seemed that you were taking $$\frac d{dx}e^{U_k}=e^{U_k}\frac{dU_k}{dx}$$ Which is wrongheaded because you really need $$\frac d{dU_j}e^{U_i}=e^{U_i}\delta_{ij}$$ Writing a simpleminded code with that change I was able to reproduce Yuriy S's analytical results.

% nonlinear.m

clear all;
close all;
m = 200;
h = 1/(m+1);
G0 = (diag(ones(m-1,1),-1)-2*diag(ones(m,1),0)+diag(ones(m-1,1),1))/h^2;
x = [1:m]/(m+1);
iarray = [1 20];
for k = 1:length(iarray),
    Uk = iarray(k)*x.*(1-x);
    err = ones(size(Uk));
    while norm(err) > 1.0e-6,
        G = G0*Uk'+exp(Uk)';
        J = G0+diag(exp(Uk));
        err = (J\G)';
        Uk = Uk-err;
    end
    U(k,:) = Uk';
end
plot([0 x 1],[0 U(1,:) 0],[0 x 1],[0 U(2,:) 0]);

Answer 2

Not related to your numerical algorithm (sorry, not my forte), but might be useful when checking the solutions. This equation is exactly solvable.

$$u''+e^u=0$$

Using inverse functions, we have for $x(u)$:

$$\frac{d^2 u}{dx^2}=-\frac{x''}{x'^3}$$

Now the equation becomes:

$$x''-e^u x'^3=0$$

New function:

$$x'=y(u)$$

Now the equation becomes:

$$y'=e^u y^3$$

Integrating:

$$\frac{1}{y^2}=C_1 - 2e^u$$

$$y = \pm \frac{1}{\sqrt{C_1-2e^u}}$$

Now returning to $x(u)$:

$$x'=\pm \frac{1}{\sqrt{C_1-2e^u}}$$

$$x=\pm \int \frac{du}{\sqrt{C_1-2e^u}}+C_2$$

The integral has an exact form:

$$\int \frac{du}{\sqrt{C_1-2e^u}}=- \frac{2}{\sqrt{C_1}} \tanh^{-1} \sqrt{1-\frac{2}{C_1} e^u}$$

So our solution is:

$$x=\pm \frac{2}{\sqrt{C_1}} \tanh^{-1} \sqrt{1-\frac{2}{C_1} e^u}+C_2$$

Or:

$$\tanh^{-1} \sqrt{1-\frac{2}{C_1} e^u}= C_2 \pm \frac{\sqrt{C_1}}{2} x$$

$$\sqrt{1-\frac{2}{C_1} e^u}=\tanh \left( C_2 \pm \frac{\sqrt{C_1}}{2} x\right)$$

$$1-\frac{2}{C_1} e^u=\tanh^2 \left( C_2 \pm \frac{\sqrt{C_1}}{2} x\right)$$

$$e^u=\frac{C_1}{2} \left(1-\tanh^2 \left( C_2 \pm \frac{\sqrt{C_1}}{2} x\right) \right)$$

$$u= \log \left( \frac{C_1}{2} \left(1-\tanh^2 \left( C_2 \pm \frac{\sqrt{C_1}}{2} x\right) \right) \right)$$

I have checked this solution by differentiating twice in wolfram alpha and it works.

Now we can try to substitute boundary conditions and get the values for the constants. From $u(0)=0$ we get:

$$C_1=\frac{2}{1- \tanh^2 C_2}$$

For the second constant we might have to solve the equation numerically. From $u(1)=0$ and the constant we just found we have:

$$1- \tanh^2 \left( C_2 \pm \frac{\sqrt{C_1}}{2}\right)=1- \tanh^2 C_2$$

Since $C_1 \neq 0$, we can use the fact that $\tanh^2$ is even to write:

$$\frac{\sqrt{C_1}}{2}-C_2=C_2$$

Or, substituting $C_1$:

$$4 C_2= \sqrt{\frac{2}{1- \tanh^2 C_2}}$$

Wolfram Alpha gives two numerical solutions:

$$C_2=2.7346756930305266999 \\ C_2 = 0.37929114976268859213$$

Now we can see how the solutions look (red the first choice for $C_2$, blue - the second choice):

Apparently, there's two solutions, both of which satisfy the boundary conditions. Which might be one of the problems with numerical algorithm.