third order nonlinear differential equation numerical solution

Hi all, I am trying to solve the following equation,

y * y''' = -1
y(0) = 1, y'(0) = 0, y''(0) = 0

to do so, I am using the following simple code with ode45

function test()

      [T, Y] = ode45(@linearized, [0 3], [1 0 0 0]);

    figure; plot(T, Y(:,1))
    axis([0 3 0 1.2])

    % linearization
    function dy = linearized(t,y)
       dy = zeros(4,1);
       dy(1) = y(2);
       dy(2) = y(3);
       dy(3) = y(4);
       dy(4) = -1*y(3)*y(1)-1;
    end
end

However, it gives me the wrong answer. Do you have any suggestion here?

Answer 1

[John D'Errico]

So lets look at how you would convert this to a system of ODEs. We have:

y*y''' = -1

So as a system, we have three variables to solve for, y, y', and y''. As you can see, you were given starting values for all three "variables", but there is no starting value for y'''. You don't need one.

I'll write out what the system is, as a system of differential equations.

y'   = y(2)
y''  = y(3)
y''' = -1/y(1)

So now we can write the function in a simple function handle form...

dy = @(t,y) [y(2);y(3);-1./y(1)];

If you prefer, an explicit nested function would have worked too, as you did, or an external m-file function.

[T, Y] = ode45(dy, [0 3], [1 0 0]);
plot(T,Y,'-')
grid on

The blue curve is y(t). It looks like something of interest happens near 1.7779 or so. If we look at the curve in blue, that is where y(t) wants to cross through zero. That would clearly be a derivative singularity, since y'''=-1/y.

[Mazi]

Thanks John for the response. Your solution is right, however when I try to use the function handle it gives me an error,

so this is the final code :

 function test()
[T, Y] = ode45(@linearized, [0 3], [1 0 0]);
plot(T,Y,'-')
grid on

    function dy = linearized(t,y)
        dy = @(t,y) [y(2);y(3);-1./y(1)];
    end
end

and it says:

Error using odearguments (line 93) TEST2/LINEARIZED returns a vector of length 1, but the length of initial conditions vector is 3.

what I am missing here?

[John D'Errico]

Because you returned a function handle there, not a vector of differential elements!

As you wrote it, do this instead:

dy = [y(2);y(3);-1./y(1)];

See that I created a function handle for my function, and then passed that directly into ode45. Whereas you have used a nested function. It needs to return a vector, which is what my solution did.