Euler method and bisection method

Question

I'd like to solve the equation $$ \phi''(x) = \lambda \sin (\phi(x)) $$ where $x \in (0,L)$, $\phi'(0) = 0$, $\phi'(L) = 0$.

Let $ \psi = \phi'$ and $$ \phi'(x) - \psi(x) = 0$$ $$ \psi'(x) - \lambda \sin (\phi(x)) = 0$$

for $x \in (0,L)$ and $\psi(0) = 0, \phi(0) = \phi_0.$

Can anybody help me to find $\phi_0 $ numerically such that $\phi'(L) = 0$ holds?

I received the advice to compute the solution of $(\phi, \psi)$ with the explicit Euler method for $\phi_0 = 1.5$ and $\phi_0 = 3$ and to use the method of bisection to compute $\phi_0$.

In addition, I received the following values:

• Number of steps (bisection): $2^6$
• Length of steps (Euler): L/100
• L = 5
• $\lambda$ = 2

Thanks for any help!

EDIT:

In the meantime, I coded a bit. I added your functions as well as an implementation of Euler and bisection. See what I did so far.

Now my problem is to connect your functions with my functions. Can you please help a bit? (For example, it's not clear to me where to define the function, and it's not clear to me when calling your functions "model" and "omegaL"...)

funtion x = eubisect()
    u = bisection(f, a, b, N, eps_step, eps_abs)

function dotu = model(t,u)
    lambda = 2;
    dotu = [ u(2); lambda*sin(u(1)) ]
end

function omegaL= f(phi0)
    L = 5;
    N = 100;
    t,u = Euler(model, 0, L, N, [phi0,0])
    omegaL = u(end,2)
end

function [t, y] = Euler(f, a, b, N, y0)
    clear t % Clears old time steps and
    clear y % y values from previous runs
    %a=0; % Initial time
    %b=1; % Final time
    %N=10; % Number of time steps
    %y0=0; % Initial value y(a)
    h=(b-a)/N; % Time step
    t(1)=a;
    y(1)=y0;
    for n=1:N % For loop, sets next t,y values
        t(n+1)=t(n)+h;
        y(n+1)=y(n)+h*f(t(n),y(n)); % Calls the function f(t,y)=dy/dt
    end
    %plot(t,y)
    %title(['Euler Method using N=',num2str(N),' steps'])
end


function [ r ] = bisection( f, a, b, N, eps_step, eps_abs )
    % Check that that neither end-point is a root
    % and if f(a) and f(b) have the same sign, throw an exception.

    if ( abs(f(a)) < eps_abs )
    r = a;
    return;
    elseif ( abs(f(b)) < eps_abs )
    r = b;
    return;
    elseif ( f(a) * f(b) > 0 )
        error( 'f(a) and f(b) do not have opposite signs' );
    end

    % We will iterate N times and if a root was not
    % found after N iterations, an exception will be thrown.

    for k = 1:N
        % Find the mid-point
        c = (a + b)/2;

        % Check if we found a root or whether or not
        % we should continue with:
        %          [a, c] if f(a) and f(c) have opposite signs, or
        %          [c, b] if f(c) and f(b) have opposite signs.

        if ( abs(f(c)) < eps_abs )
            r = c;
            return;
        elseif ( f(c)*f(a) < 0 )
            b = c;
        else
            a = c;
        end

        % If |b - a| < eps_step, check whether or not
        %       |f(a)| < |f(b)| and |f(a)| < eps_abs and return 'a', or
        %       |f(b)| < eps_abs and return 'b'.

        if ( b - a < eps_step )
            if ( abs( f(a) ) < abs( f(b) ) && abs( f(a) ) < eps_abs )
                r = a;
                return;
            elseif ( abs( f(b) ) < eps_abs )
                r = b;
                return;
            end
        end
    end

    error( 'the method did not converge' );
end

Answer 1

You are to apply the single shooting method.

You already got the first order system

function dotu = model(t,u)
    lambda = 2;
    dotu = [ u(2); lambda*sin(u(1)) ]
end

Now define the shooting function

function omegaL= f(phi0)
    L = 5;
    t,u = ode45(model, [0,L], [phi0,0])
    omegaL = u(end,2)
end

and then call your bisection routine,

phi0 = bisection(f, 1.5, 3.0, 64).

---

Graph of function $f$, the root inside the given interval is close to $\phi_0=1.80132467093$. This is for a fairly exact integration of the IVP. For the rather inexact Euler method you will get a distorted graph and thus different roots.

Contents of bvp_pendulum, this is for octave, replace endif, endfor, endfunction with end if necessary.

function phi0 = bvp_pendulum()
    phi0 = bisection(@f, 1.5, 3.0, 64, 1e-10, 1e-10)
endfunction


function dotu = model(t,u)
    lambda = 2;
    dotu = [ u(2), lambda*sin(u(1)) ];
endfunction;

function [t,y] = Euler(f, a, b, N, y0)
    h=(b-a)/N; % Time step
    t(1)=a; % Initial time
    y(1,:)=y0; % Initial value y(a)
    for n=1:N % For loop, sets next t,y values
        t(n+1)=t(n)+h;
        y(n+1,:)=y(n,:)+h*f(t(n),y(n,:)); % Calls the function f(t,y)=dy/dt
    endfor;
endfunction;

function omegaL = f(phi0)
    L = 5; N=100;
    [t,u] = Euler(@model, 0, L, N, [phi0,0]);
    omegaL = u(end,2);
endfunction;

function r = bisection( f, a, b, N)
    % if f(a) and f(b) have the same sign, throw an exception.
    fa = f(a);
    fb = f(b)
    if ( fa * fb > 0 )
        error( 'f(a) and f(b) do not have opposite signs' );
    endif

    % We will iterate N times

    for k = 1:N
        % Find the mid-point
        c = (a + b)/2;

        % we continue with:
        %          [a, c] if f(a) and f(c) have opposite signs, or
        %          [c, b] if f(c) and f(b) have opposite signs.
        fc = f(c);
        disp([int2str(k),': new point c=',num2str(c,'%20.17f'),', f(c)=',num2str(fc,'%20.17f')]);
        if ( fc*fa < 0 )
            b = c; fb = fc;
        else
            a = c; fa = fc;
        endif

    endfor
    r = c;
endfunction