3
$\begingroup$

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
$\endgroup$
  • $\begingroup$ This is the physical pendulum seen from the upper, unstable equilibrium point. The stable point is at $ϕ=\pi$ so that for $ϕ_0\approx\pi$ you get period $2\pi$ and for $ϕ_0\approx0$ you get arbitrarily large periods. $L$ has to be a multiple of the half-period, so $L=5>\pi$ is admissible. $\endgroup$ – LutzL May 16 '18 at 12:24
  • $\begingroup$ Yeah very cool! You knew this, but how is it possible to see the answer without knowing it in advance? Just try some values? (If yes, do you have any code on Matlab ?) $\endgroup$ – DMan May 16 '18 at 12:37
  • $\begingroup$ These are basic facts on the pendulum equation. What do you know about the single shooting method? In matlab you would normally use the boundary value problem solver, but then there would be no bisection method be needed. $\endgroup$ – LutzL May 16 '18 at 12:53
  • $\begingroup$ Ok :) No, by "using Matlab", I meant that I implement the Euler and bisection methods and that I try, try and try to get $\phi_0 = \pi$ :) $\endgroup$ – DMan May 16 '18 at 13:17
  • $\begingroup$ Then please document what you already programmed and computed. $\endgroup$ – LutzL May 16 '18 at 13:54
0
$\begingroup$

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. enter image description here

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
$\endgroup$
  • $\begingroup$ Thanks for this answer! I added an edit, maybe you can check this out? Thanks a lot! $\endgroup$ – DMan May 17 '18 at 7:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.