Trying to solve a two-point boundary value problem on MATLAB

Question

I have the following ODE

$$ u'' = -(1 + e^u), u(0)=0, u(1)=1$$

with boundary conditions and $t \in (0,1)$. I want to solve this ODE using the shooting method. First,we convert this to a system of equations. Let $y_1=u, y_2=u'$ and so

$$ \begin{bmatrix} y_1' \\ y_2' \end{bmatrix} = \begin{bmatrix} y_2 \\ -(1+e^{y_1}) \end{bmatrix}, \; \; \; \; \begin{bmatrix} y_1(0) \\ y_2(0) \end{bmatrix} = \begin{bmatrix} 0 \\ u'(0) \end{bmatrix}$$

Im trying to understand how the shooting method works and hopefully code it in matlab. So,I have to make a guess for $u'(0)$ and then use a ODE solve to get $u(1)=1$. Is this how we procced with this method?

Answer 1

As you did, you first write the system as

$\begin{bmatrix} y_1' \\ y_2' \end{bmatrix} = \begin{bmatrix} y_2 \\ -(1+e^{y_1}) \end{bmatrix}, \; \; \; \; \begin{bmatrix} y_1(0) \\ y_2(0) \end{bmatrix} = \begin{bmatrix} 0 \\ u'(0) \end{bmatrix}$

in order to use a standard numerical method for ODEs.

So, the goal is to give the right value to $u'(0)$ in order to have $u(1)=1$. In your system, you have to guess the value of $y_2(0)$.

Of course you could do it by hand, but you can actually write a bisection method in order to find the correct value for the first derivative. The function you have to set to zero is, actually, $F=y_1(N)-1$.

So, assume to use a numerical method for solving your ODE, say explicit Euler (EE) (for sake of simplicity):

• Consider two initial data $y_0= [0,a]^T$, $y_0=[0,b]^T$
• Solve the ODE with (EE) those different initial values. Then you compute, for both the numerical solutions you get $f_a=y_1(N)-1$ and $f_b=y_1(N)-1$.
• Compute the mean of $a,b$, call it $x_m$, and set another initial data $y_0=[0,x_m]^T$, and solve the ODE with (EE) and the last initial data and define $f_m=y_1(N)-1$.
• Apply a standard bisection method: if $f_m \cdot f_a <0$, then $b=x_m$ and $f_b=f_m$, otherwise $a=x_m,f_a=f_m$
• Repet the above procedure until the distance between $a$ and $b$ is $|b-a| \leq \epsilon$, where $\epsilon$ is a chosen tolerance

The following (runable) MatLab code computes the right value for $y_2(0)$, which is approximately $2.4592$.

fun=@(t,y) [y(2);-1-exp(y(1))];  
m=200; %time steps
tstar=1;
k=tstar/m;
t = linspace(0,tstar,m+1);
y=zeros(2,m+1);
a=-5;
b=5; % extreme initial points for bisection
iter=0;
maxit=100;
tol=1e-6; %tolerance

while (abs(b-a)>tol) && (iter<maxit)

y(:,1)=[0,a]; %a initial guess iniziale for y2

for n = 1:m
    y(:,n+1) = y(:,n)+k*fun(t(n),y(:,n)); %EE
end
fa=y(1,end)-1; %I have to set to zero F=y(1) - 1. (Want y(1)=1)

y(:,1)=[0,b];

for n = 1:m
    y(:,n+1) = y(:,n)+k*fun(t(n),y(:,n)); %EE
end
fb=y(1,end)-1;    


iter=iter+1;
xm=(a+b)/2;
y(:,1)=[0,xm];

for n = 1:m
    y(:,n+1) = y(:,n)+k*fun(t(n),y(:,n)); %EE
end

fm=y(1,end)-1;


if fa*fm<=0
    b=xm;
    fb=fm;
else 
    a=xm;
    fa=fm;
end

end


disp('The correct value is')
disp(xm)

%%%PLOT THE NUM SOLUTION%%

y=zeros(2,m+1);
y(:,1)=[0,2.4592];
for n = 1:m
     y(:,n+1) = y(:,n)+k*fun(t(n),y(:,n)); %EE
end
plot(t,y(1,:),'-*')

The correct numerical solution is the following