solve a system of ode by bvp4c

Hi, I am trying so solve a system of boundary value problem. I have tried looking into the examples given yet I still have problem in writing the code and running the program.

Briefly, I have 3 equations to solve. I have convert them to first order system in order to put it in matlab. Where, I let W=y1, W'=y2, U=y3, U'=y4, U''=y5, V=y6, V'=y7, V''=y8. I have the boundary 0f x from [0 infinity], and I choose infinity=20. I also need to plot graphs of x-axis:x (from 0 to 20) for y-axis:U, V and W respectively.

I know there are lots of mistakes here hence please correct me if Im wrong. I have attached the code here.

Thank you :))

[Torsten]

Please include equations and boundary conditions in a mathematical notation.

Best wishes

Torsten.

[Day Rosli]

So this is the equation and boundary conditions:

2*U + x*(1-n/n+1)*U' + W' = 0

U^2 - (V+1)^2 + [W+(1-n/n+1)*x*U]*U' - {(U'^2+V'^2)^((n-1)/2)*U''+U'*[(n-1)*(U'U''+V'V'')*(U'^2+V'^2)^((n-3)/2)]}=0

2*U*(V+1) + [W+(1-n/n+1)*x*U]*V' - {(U'^2+V'^2)^((n-1)/2)*V''+V'*[(n-1)*(U'U''+V'V'')*(U'^2+V'^2)^((n-3)/2)]}=0

bc:

U(0)=V(0)=W(0)=0, U(infinity)=0, V(infinity)=-1.

So I transformed them into 1st order by letting

W=y1, W'=y2, U=y3, U'=y4, U''=y5, V=y6, V'=y7, V''=y8.

Thus, I have this equation.

y(2) = -2*y(3)-x*(1-n/n-1)*y(4);

y(4) = (-y(3)^2+(y(6)+1)^2+((y(4)^2+y(7)^2)^((n-1)/2)*y(5)+y(4)*((n-1)*(y(4)*y(5)+y(7)*y(8))*(y(4)^2+y(7)^2)^((n-3)/2))))/(y(1)+(1-n/n+1)*x*y(3))

y(7) = (-2*y(3)*(y(6)+1)+((y(4)^2+y(7)^2)^((n-1)/2)*y(8)+y(7)*((n-1)*(y(4)*y(5)+y(7)*y(8))*(y(4)^2+y(7)^2)^((n-3)/2))))/(y(1)+(1-n/n+1)*x*y(3))

bc: f(1)(0)=f(3)(0)=f(6)(0)=0, f(3)(infinity)=0, f(6)(infinity)=-1.

[Taylor Nichols]

how did you arange the boundary conditions in your code if you don't mind me asking?

I've only done a ODE with one dependent variable in which I arranged my conditions this way

This is a 4th order equation btw.

bc = [yl(1) - hi; 
      yl(2);
      yr(1) - ho;
      yr(2)];
    
        

Answer 1

[Torsten]

Solve your equations for W', U'' and V''.

Let

y(1) = W, y(2) = U, Y(3) = U', Y(4) = V, Y(5) = V'.

Your system then reads

Y(1)' = F0(Y(1),Y(2),Y(3),Y(4),Y(5))
Y(2)' = Y(3)
Y(3)' = F1(Y(1),Y(2),Y(3),Y(4),Y(5))
Y(4)' = Y(5)
Y(5)' = F2(Y(1),Y(2),Y(3),Y(4),Y(5))

Now use one of the ODE solvers in the usual manner.

Best wishes

Torsten.

[Day Rosli]

Thank you so much for this. However, I have problem to proceed for Y(3)' since there is V'' in the equations. Same goes to Y(5)'. So, do you have any idea about this?

Let

Y(1) = W, Y(2) = U, Y(3) = U', Y(4) = V, Y(5) = V'

Hence,

 Y(1)' = -2Y(2)+x(1-n/n+1)Y(3)
Y(2)' = Y(3)
Y(3)' = Y(2)^2-(Y(4)+1)^2+(Y(1)+x(1-n/n+1)Y(2))Y(3)-Y(3)(n-1)(Y(3)^2+Y(5)^2)^((n-3)/2)Y(5)V''/(Y(3)^2+Y(5)^2)^((n-1)/2)+Y(3)(n-1)(Y(3)^2+Y(5)^2)((n-3)/2)
Y(4)' = Y(5)
Y(5)' = -2Y(2)(Y(4)+1)+(Y(1)+x(1-n/n+1)Y(2))Y(5)-Y(5)(n-1)(Y(3)^2+Y(5)^2)^((n-3)/2)Y(3)U''/(Y(3)^2+Y(5)^2)^((n-1)/2)+Y(5)(n-1)(Y(3)^2+Y(5)^2)((n-3)/2)

So, as you can see as above, at Y(3)', there exist V'' and also at Y(5)', there exist U''. So how do I assign that?

Thank you in advance.

[Torsten]

Order equations (2) and (3) such that they look like

a1*U'' + b1*V'' = c1
a2*U'' + b2*V'' = c2

with a1, b1, c1, a2, b2, c2 expressions only depending on U, U', V, V', W and W'.

Then solve this (linear) system of equations from above for U'' and V'':

U'' = (c1*a2-c2*a1)/(b1*a2-b2*a1)
V'' = -(c1*b2-c2*b1)/(b1*a2-b2*a1)

Best wishes

Torsten.

P.S.: It seems you forget parentheses around expressions, e.g. (1-n/n+1) should be (1-n/(n+1)), I guess.

[Day Rosli]

Thank you again, for your help. Unfortunately, I am not very clear of some of the explanation above (Sorry, my bad). Okay, from what I understand so far, I have to write the equations so that it will look like this, right?

a1*U'' + b1*V'' = c1
a2*U'' + b2*V'' = c2

Okay, so this is what I got from rearranging equations for Y(3)' and Y(5)'

[((U'^2+V'2)^((n-1)/2)+U'(n-1)(U'^2+V'^2)^((n-3)/2)U']U''+[U'(n-1)(U'^2+V'^2)^((n-3)/2)V']V''=U^2-(V+1)^2+(W+x(1-n/(n+1))U)U'

[V'(n-1)(U'^2+V'^2)^((n-3)/2)U']U''+[((U'^2+V'2)^((n-1)/2)+V'(n-1)(U'^2+V'^2)^((n-3)/2)V']V''=-2U(V+1)+(W+x(1-n/(n+1))U)V'

in which

a1=((U'^2+V'2)^((n-1)/2)+U'(n-1)(U'^2+V'^2)^((n-3)/2)U' 
b1=U'(n-1)(U'^2+V'^2)^((n-3)/2)V'
c1=U^2-(V+1)^2+(W+x(1-n/(n+1))U)U'
a2=V'(n-1)(U'^2+V'^2)^((n-3)/2)U'
b2=((U'^2+V'2)^((n-1)/2)+V'(n-1)(U'^2+V'^2)^((n-3)/2)V'
c2=-2U(V+1)+(W+x(1-n/(n+1))U)V'

so, what do you mean by

'Then solve this (linear) system of equations from above for U'' and V'':
U'' = (c1*a2-c2*a1)/(b1*a2-b2*a1)
V'' = -(c1*b2-c2*b1)/(b1*a2-b2*a1)'  ?

Kindly need your help, thank you in advance.

Answer 2

[Torsten]

I mean that the function F1 from above is given by

F1(W,U,U',V,V') = (c1*a2-c2*a1)/(b1*a2-b2*a1)

and the function F2 from above by

F2(W,U,U',V,V') = -(c1*b2-c2*b1)/(b1*a2-b2*a1)

Now you can define your array "dydx" in function "ex11ode" by

dydx = [-2Y(2)+x(1-n/n+1)Y(3)
        Y(3)
        F1(Y(1),Y(2),Y(3),Y(4),Y(5))
        Y(5)
        F2(Y(1),Y(2),Y(3),Y(4),Y(5))]

Best wishes

Torsten.

[Day Rosli]

Okay, now I got what u meant and then when I tried to run the program, it gives the new error which is

Unable to solve the collocation equations -- a singular Jacobian encountered.

So, what does this means? Does this has something to do with my equation? And do you mind to explain what should I do?

Thank you so much.

function ex111bvp
solinit = bvpinit(linspace(0,20,2),[0 0 0 0 -1]);
sol = bvp4c(@ex111ode,@ex111bc,solinit);
solinit=sol;
figure
plot(sol.x,sol.y(2,:));
title('Example 111')
ylabel('U')
xlabel('x')
% --------------------------------------------------------------------------
function dydx = ex111ode(x,y,n)
n=1;
dydx =  [ -2*y(3)-x*(1-n/(n+1))*y(4)
        y(3)
        y(3)*y(5)*(y(3)^2+y(5)^2)^((n-3)/2)*(n-1)*(y(3)*(y(1)-x*y(2)*(1-n/(n+1)))-(y(4)+1)^2+y(2))-((y(3)^2+y(5)^2)^((n-1)/2)+y(3)^2*(y(3)^2+y(5)^2)^((n-3)/2)*(n-1))*((y(3)^2+y(5)^2)^((n-1)/2)+y(5)*(y(3)^2+y(5)^2)^((n-1)/2)+y(5)*(y(3)^2+y(5)^2)^((n-3)/2)*(n-1))-((y(3)^2+y(5)^2)^(3-n)*(y(5)*(y(1)*x*y(2)*(1-n/(n+1)))-2*y(2)*(y(4)+1))*((y(3)^2+y(5)^2)^((n-1)/2)+y(3)*(y(3)^2+y(5)^2)^((n-3)/2)*(n-1))) 
        y(5)
        ((y(3)^2+y(5)^2)^((n-1)/2)+y(3)^2*(y(3)^2+y(5)^2)^((n-3)/2)*(n-1))*((y(3)^2+y(5)^2)^((n-1)/2)+y(5)^2*(y(3)^2+y(5)^2)^((n-3)/2)*(n-1))-((y(3)^2+y(5)^2)^((n-1)/2)+y(5)^2*(y(3)^2+y(5)^2)^((n-3)/2)*(n-1))*(y(3)*(y(1)-x*y(2)*(1-n/(n+1)))-(y(4)+1)^2+y(2))+((y(3)^2+y(5)^2)^(3-n)*(y(3)^2+y(5)^2)^((n-3)/2)*(y(5)*(y(1)-x*y(2)*(1-n/(n+1)))-2*y(2)*(y(4)+1)))/(y(3)*y(5)*(n-1))];
%-------------------------------------------------------------------------
function res = ex111bc(ya,yb)
res = [ ya(1)
      ya(2)
      ya(4)
      yb(2)
      yb(4)+1];
%-------------------------------------------------------------------------

[Torsten]

Answered under

http://www.mathworks.com/matlabcentral/answers/264072-error-using-bvp4c-singular-jacobian-encountered

Best wishes

Torsten.

[Day Rosli]

Hi, Torsten. I manage to run this code but up until infinity=10. However, I have to plot until infinity=20. At first, when I tried infinity=10, it works well but when I wanted to extend infinity=20, the graphs looks really bad. The solution should be approaches to 0 as x=20 (asymptotic).

Hence, I looked at some example on continuation where I ended at infinity=10 and use continuation for infinity>10 to 20. However, it did not work.

This is without continuation:

function Script
n=1;
infinity=10;
maxinfinity=20;

solinit = bvpinit(linspace(0,infinity,100),[0 0 0 0 -1]);
sol = bvp4c(@(x,y)VK(x,y,n),@VKbc,solinit);
figure
  plot(sol.x,sol.y(3,:),sol.x(end),sol.y(3,:))
  title('Example Script')
  axis([0 maxinfinity 0 -1.5]);
  ylabel('U')
  xlabel('eta')
  hold on
  drawnow
  shg

function yprime = VK(x,y,n)
n=1;
a = ((1-n)/(n+1))*x;
c = (y(4)^2+y(5)^2)^((1-n)/(n+1));
yprime = [ c*y(4)
           c*y(5)
           -2*y(1) - a*c*y(4)
           y(1)^2 - (y(2)+1)^2 + (y(3)+a*y(1))*c*y(4)
           2*y(1)*(y(2)+1) + (y(3)+a*y(1))*c*y(5)];
end

function res = VKbc(ya,yb)

res = [ya(1);ya(2);ya(3);yb(1);yb(2)+1];

end

end

And this is when I use continuation where nothing happen:

function Script
infinity=10;
maxinfinity=20;
solinit = bvpinit(linspace(0,infinity,800),[0 0 0 0 -1]);
sol = bvp4c(@VK,@VKbc,solinit);
figure
  plot(sol.x,sol.y(3,:),sol.x(end),sol.y(3,:))
  title('Example Script')
  axis([0 maxinfinity 0 -1.5]);
  ylabel('U')
  xlabel('eta')
  hold on
  drawnow
  shg
    for Bnew = infinity+1:maxinfinity
     solinit = bvpxtend(sol,Bnew);   % Extend the solution to Bnew.
     sol = bvp4c(@VK,@VKbc,solinit);
     plot(sol.x,sol.y(3,:),sol.x(end),sol.y(3,:))
     drawnow
    end
  hold off
Script
function yprime = VK(x,y,n)
n=1;
a = ((1-n)/(n+1))*x;
c = (y(4)^2+y(5)^2)^((1-n)/(n+1));
yprime = [ c*y(4)
           c*y(5)
           -2*y(1) - a*c*y(4)
           y(1)^2 - (y(2)+1)^2 + (y(3)+a*y(1))*c*y(4)
           2*y(1)*(y(2)+1) + (y(3)+a*y(1))*c*y(5)];
end
function res = VKbc(ya,yb)
res = [ya(1);ya(2);ya(3);yb(1);yb(2)+1];
end
end

I referred to this example . So, if you have time, mind to help me out? Thank you in advance!

[Torsten]

Why don't you just solve the ODEs with right limit x=20 instead of x=10 ?

Does the solver encounter difficulties finding a solution ?

Best wishes

Torsten.

[Day Rosli]

Hi Torsten. I have made a mistake regarding the latest question I asked. there is no problem when running the code above. However, I also need to plot c against x. Where c is a function of

c = (y(4)^2+y(5)^2)^((1-n)/(n+1))

as stated on the coding above. Do you know how?

[Torsten]

Since you set n=1, c is identically 1.

Best wishes

Torsten.

[Day Rosli]

Yes, I know that. But what if n is other values? I mean how should I write a code to insert in the main code in order to plot the function c?

[Torsten]

Something like

c=(sol.y(4,:).^2+sol.y(5,:).^2).^((1-n)/(n+1));
plot(sol.x,c);

after the call to bvp4c should work.

Best wishes

Torsten.

[Day Rosli]

Thank you Torsten! It works well.