MATLAB Answers

How to solve following BVP ODE?

3 views (last 30 days)
Hi mates,
I am trying to solve the following ODE in MATLAB, but I do not attain a reasonable answer. It would be highly appreciated if you let me know of the bugs in my code below the ODE.
The ODE is: (d/dx)(y^3 dy/dx)+(2/3)(x dy/dx)-(1/3)y=0;
The BC's are: (y^3)(dy/dx)=-1, for x=0; and y(inf)=0;
What I have done is to use transforms as below
Y(1)=y;
Y(2)=dy/dx;
Then, dY(1)/dx=Y(2) and dY(2)/dx=((1/3)(1/Y(1)^2))-((2/3)(xY(2)/Y(1)^3)) (to have this equation I just ignored the high order small value of (dy/dx)^2 which is an assumption I took; If there is better solution I would be happy to replace).
According to above descriptions, I used following code in MATLAB:
function liftoff()
xmesh=linspace(0,4,10);
solinit = bvpinit(xmesh,[1 0]);
sol = bvp4c(@lode,@bcs,solinit);
plot(sol.x,sol.y(1,:),'b-x');
function dYdx = lode(x,Y)
dYdx(1) = Y(2);
dYdx(2) = (1/3)*((1/(Y(1)^2))-(2*x*Y(2)/(Y(1)^3)));
function res = bcs(ya,yb)
res = [ ((ya(1)^3)*ya(2))+1;
yb(1)];
But, I encountered with error of "a singular Jacobian".
I am looking forward to hearing your advices.
Cheers,
Dariush

Accepted Answer

Alan Stevens
Alan Stevens on 13 May 2021
Weird ODE and initial boundary condition! You can solve for small values of x using the following (where y(x=0) was chosen by trying a few values). Larger values of y0 lead to an upturn in y after the dive!
y0 = 1.311435;
v0 = -1./y0.^3;
Y0 = [y0, v0];
xend = 4;
xspan = [0 xend];
[x, Y] = ode15s(@fn, xspan, Y0);
y = Y(:,1);
v = Y(:,2);
plot(x,y),grid
function dYdx = fn(x,Y)
y = Y(1);
v = Y(2);
if y<=0 % just to prevent weird behaviour !!
dYdx = [0; 0];
else
dYdx = [v;
(y/3 - 2*x.*v/3 - 3*y.^2.*v.^2 )./y.^3];
end
end
  1 Comment
Dariush Javani
Dariush Javani on 13 May 2021
Thanks Alan for the answer! It is exactly what I was trying to do.
cheers,
Dariush

Sign in to comment.

More Answers (0)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!