Solving a third order non-linear ode using ode45

I need to solve

F''' + 2FF" + 1 - F'^2 = 0

with the boundary conditions

F(0)=F'(0)=0 and F'(infinity)=1.

I am new to using the ode solver in matlab and am not sure how to make it solve a non-linear third order equation. Any suggestion would be appreciated.

[Shiv Patel]

Thank you for this comment, it helped with my homework!

Answer 1

[Star Strider]

The ode45 function is not applicable to two-point boundary value problems such as yours. You need to use bvp4c for example, as I did here:

dYdX = @(X,Y) [Y(2);  Y(3);  Y(2).^2-Y(1).*Y(2).*2.0-1.0];  % Differential equation

res = @(ya,yb) [ya(1);  ya(2);  yb(2)-1];                   % Boundary conditions

SolYinit = bvpinit([0 1E+1], [1; 1; 1]);

Fsol = bvp4c(dYdX, res, SolYinit);

X = Fsol.x;
F = Fsol.y;

figure(1)
plot(X, F)
legend('F_1', 'F_2', 'F_3', 3) 
grid

You will not be able to set your endpoint to Inf or other number larger than about 10 (at least I could not) because that generates a singular Jacobian (or so MATLAB tells me). There are links to bvpinit and the other functions throughout and at the end of the bvp4c page. You might also want to experiment with bvp5c.

This was an interesting problem. I am glad you asked the question.

[Abhinav]

Thanks for taking the time to go through my question. The plots being generated by thus code do not match the results that I was expecting to get, and I am not able to figure out why, because everything in the code seems good.

I finally did solve this with ODE45 and here is the code :

 clc
 clear all

 xmin=0;
 xmax=2;
 while((xmax-xmin)>10e-6)
    x=(xmin+xmax)/2;
    init = [0 0 x];
    eta = [0 5];
    [t,y] = ode45('eqn',eta,init)
    %max(y(:,2))
     if max(y(:,2))>1
         xmax=x;
     else
         xmin=x;
     end
 end

plot(t,y(:,3),'r+-',t,y(:,2),'b*-',t,y(:,1),'g+-')
axis equal
axis([0 3 0 2])
legend('d2F','dF','F')
grid on

The function was defined in another file as so....

    function f = eqn(t,y)

    f = zeros(size(y));

    f(1) = y(2);
    f(2) = y(3);
    f(3) = -1 -2*y(1)*y(3) + y(2)*y(2);

The plots that I am getting from this are exactly what I wanted.

It mainly involved figuring the initial conditions for the given boundary condition, i.e., finding a value of F"(0) so that F'(infinity)=1, which this piece of code calculates as 1.3119 approximately. Here, I tried different values for 'infinity' and finally settled on 5 as a good number.

[Star Strider]

I'm glad I could help.

In the bvp4c version, changing the bvpinit call to:

SolYinit = bvpinit([0 5], [0; 0; 0]);

produces substantially the same curves as your code, with slightly greater magnitudes by a factor of about 1.4 between 0 and 3. Setting the third initial condition to 1.3 doesn't appear to make a significant difference. I didn't see a difference in your and my differential equation functions, other than I set mine up as an anonymous function.

[Ngoc-Tram Le]

Thanks a lot for your sharing this code. I got the same problem like that.

[Saeid]

I have to solve a very similar problem:

F''' + 1 - F'^2 = 0

But my BCs are different:

F'(0)=F"(inf)=0

F'(inf)=1

I formulate the problem exactly as suggested above:

dYdX = @(X,Y) [Y(2);  Y(3);  Y(2).^2-1];  % Differential equation
res = @(ya,yb) [ya(2);  yb(3);  yb(2)-1];   

But when I try Star Strider's method it generates the Jacobian error. Could it be that the BCs do not contain F(...)?

[David After]

How can i solve this eqaution?

F''' (eta)+ F(eta)F"(eta) + F'^2(eta) = 0

with the boundary conditions

F(0)=F''(0)=0 and F'(infinity)=0 and eta is 0:0.1:6

i want to plot it and create a table for it (eta-f-f'-f'')

I am new to using the ode solver in matlab and am not sure how to make it solve a equation. Any suggestion would be appreciated.

thank you

my email : [e-mail]

[David After]

How can i solve this eqaution?

F''' (eta)+ F(eta)F"(eta) + F'^2(eta) = 0

with the boundary conditions

F(0)=F''(0)=0 and F'(infinity)=0 and eta is 0:0.1:6

i want to plot it and create a table for it (eta-f-f'-f'')

I am new to using the ode solver in matlab and am not sure how to make it solve a equation. Any suggestion would be appreciated.

thank you

my email : [e-mail]

Answer 2

[R Vishnu]

Hi I hope to correct a minor error. it is f(3) = -1 -2*y(1)*y(2) + y(2)*y(2); and not f(3) = -1 -2*y(1)*y(3) + y(2)*y(2); as per your equation F''' + 2FF' + 1 - (F'^2) = 0.

[Abhinav]

Hi Vishnu, My bad.. the second term was supposed to be 2FF". I have edited the original question too.

Thanks for pointing that out!

Answer 3

[Hemanth kumar C]

I need to solve

(D^2-a^2)fi+a*[R*thetha+(Ra*Phi/Le)] = 0,

(D^2-a^2)thetha-(a*(dt/dy)*fi)-Pe*D(thetha)+gamma*Phi=0,

(D^2-a^2)Phi-Le*Pe*D(Phi)-a*le*(dc/dy)*fi=0

with the boundary conditions y=0;fi=Phi=thetha=0

I am new to using the ode solver in matlab and am not sure how to make it solve a non-linear second order three equation and i have written the program but i am not getting proper output . if you want to see my practise paper i will attach here . Any suggestion would be appreciated.

[Torsten]

I don't understand your equations.

If fi, Phi and thetha are your dependent variables and y is the independent variables, what does dc/dy, dt/dy mean in your equations ?

Please clarify which variables are given constants, dependent and independent variables.

The best will be to delete your question (answer ?) here and open a new thread.

Best wishes

Torsten.

[Hemanth kumar C]

dc/dy ,dt/dy are the temperature and concentration gradients.i derived equations for that. i attached the program sir you can see the both equations

[Hemanth kumar C]

i am solving this above attached paper .i derived everything only programming part i am not getting few things .please help me to solve this problem. you can see what is dc/dy and dt/dy in the paper.

Answer 4

[David After]

How can i solve this eqaution?

F''' (eta)+ F(eta)F"(eta) + F'^2(eta) = 0

with the boundary conditions

F(0)=F''(0)=0 and F'(infinity)=0 and eta is 0:0.1:6

i want to plot it and create a table for it (eta-f-f'-f'')

I am new to using the ode solver in matlab and am not sure how to make it solve a equation. Any suggestion would be appreciated.

thank you

my email : [e-mail]