Solving 2nd order differential equation with boundary condition z(inf) = 0

Question

How can I solve a 2^nd order differential equation with boundary condition like z(inf)?

2(x+0.1)·z'' + 2.355·z' - 0.71·z = 0
z(0) = 1
z(inf) = 0
z'(0) = -4.805

I can't understand where the boundary value z(inf) is to be used in ode45() function. I used in the following condition [z(0) z'(0) z(inf)], but this does not give accurate output.

function [T, Y]=test()

    % some random x function 
    x = @(t) t;

    t=[0 :.01 :7];


    % integrate numerically
    [T, Y] = ode45(@linearized, t, [1 -4.805 0]);

    % plot the result
    plot(T, Y(:,1))

    % linearized ode
    function dy = linearized(t,y)

       dy = zeros(3,1);
       dy(1) = y(2);
       dy(2) = y(3);
       dy(3) = (-2.355*y(2)+0.71*y(1))/((2*x(t))+0.2);

    end

end

please help me to solve this differential equation.

Answer 1

You seem to have a fairly advanced problem on your hands, but very limited knowledge of MATLAB and/or ODE theory. I'm happy to explain more if you want, but that should be in chat (I'll invite you) or via personal e-mail (my last name AT the most popular mail service from Google DOT com)

Now that you've clarified a few things and explained the whole problem, things are a bit more clear and I was able to come up with a reasonable solution. I think the following is at least in the general direction of what you'd need to do:

function  [tSpan, Y2, Y3] = test

    %%# Parameters

    %# Time parameters
    tMax  = 1e3;
    tSpan = 0 : 0.01 : 7;

    %# Initial values
    y02 = [1 -4.805];    %# second-order ODE
    y03 = [0 0 4.8403];  %# third-order ODE

    %# Optimization options
    opts = optimset(...
        'display', 'off',...
        'TolFun' , 1e-5,...
        'TolX'   , 1e-5);


    %%# Main procedure

    %# Find X so that z2(t,X) -> 0 for t -> inf
    sol2 = fminsearch(@obj2, 0.9879680932400429, opts);

    %# Plug this solution into the original
    %# NOTE: we need dense output, which is done via deval()
    Z = ode45(@(t,y) linearized2(t,y,sol2), [0 tMax], y02);

    %# plot the result
    Y2 = deval(Z,tSpan,1);
    plot(tSpan, Y2, 'b');


    %# Find X so that z3(t,X) -> 1 for t -> inf
    sol3 = fminsearch(@obj3, 1.215435887288112, opts);

    %# Plug this solution into the original
    [~, Y3] = ode45(@(t,y) linearized3(t,y,sol3), tSpan, y03);

    %# plot the result
    hold on, plot(tSpan, Y3(:,1), 'r');

    %# Finish plots
    legend('Second order ODE', 'Third order ODE')
    xlabel('T [s]')
    ylabel('Function value [-]');

    %%# Helper functions

    %# Function to optimize X for the second-order ODE
    function val = obj2(X)
        [~, y] = ode45(@(t,y) linearized2(t,y,X), [0 tMax], y02);
        val = abs(y(end,1));
    end

    %# linearized second-order ODE with parameter X
    function dy = linearized2(t,y,X)
        dy = [
            y(2)
            (-2.355*y(2) + 0.71*y(1))/2/(X*t + 0.1)
            ];
    end


    %# Function to optimize X for the third-order ODE
    function val = obj3(X3)
        [~, y] = ode45(@(t,y) linearized3(t,y,X3), [0 tMax], y03);
        val = abs(y(end,2) - 1);
    end

    %# linearized third-order ODE with parameters X and Z
    function dy = linearized3(t,y,X)
        zt = deval(Z, t, 1);
        dy = [
            y(2)
            y(3)
            (-1 -0.1*zt + y(2) -2.5*y(3))/2/(X*t + 0.1)
            ];
    end

end

Answer 2

As in my comment above, I think you're confusing a couple of things. I suspect this is what is requested:

function [T,Y] = test

    tMax = 1e3;

    function val = obj(X)
        [~, y] = ode45(@(t,y) linearized(t,y,X), [0 tMax], [1 -4.805]);
        val = abs(y(end,1));
    end

    % linearized ode with parameter X
    function dy = linearized(t,y,X)
        dy = [
            y(2)
            (-2.355*y(2) + 0.71*y(1))/2/(X*t + 0.1)
            ];
    end

    % Find X so that z(t,X) -> 0 for t -> inf
    sol = fminsearch(@obj, 0.9879);

    % Plug this ssolution into the original 
    [T, Y] = ode45(@(t,y) linearized(t,y,sol), [0 tMax], [1 -4.805]);

    % plot the result
    plot(T, Y(:,1));

end

but I can't get it to converge anymore for tMax beyond 1000 seconds. You may be running into the limits of ode45 capabilities w.r.t. accuracy (switch to ode113), or your initial value is not accurate enough, etc.