solve a system of ode by bvp4c

Question

Day Rosli 2016-1-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/263650-solve-a-system-of-ode-by-bvp4c

评论： Taylor Nichols 2019-1-24

ex111bvp.m

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 :))

3 个评论
显示 1更早的评论隐藏 1更早的评论

Day Rosli 2016-1-13

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 2019-1-24

在 MATLAB Online 中打开

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 2016-1-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/263650-solve-a-system-of-ode-by-bvp4c#answer_206249

在 MATLAB Online 中打开

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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Torsten 2016-1-15

编辑：Torsten 2016-1-15

在 MATLAB Online 中打开

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 2016-1-16

在 MATLAB Online 中打开

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 2016-1-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/263650-solve-a-system-of-ode-by-bvp4c#answer_206493

在 MATLAB Online 中打开

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.

9 个评论
显示 7更早的评论隐藏 7更早的评论

Day Rosli 2016-1-18

在 MATLAB Online 中打开

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];
%-------------------------------------------------------------------------

Day Rosli 2016-1-19

编辑：Day Rosli 2016-1-19

在 MATLAB Online 中打开

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 2016-1-25

在 MATLAB Online 中打开

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 2016-1-25

Thank you Torsten! It works well.

请先登录，再进行评论。

solve a system of ode by bvp4c

3 个评论
显示 1更早的评论隐藏 1更早的评论

采纳的回答

3 个评论
显示 1更早的评论隐藏 1更早的评论

更多回答（1 个）

9 个评论
显示 7更早的评论隐藏 7更早的评论

另请参阅

类别

标签

Community Treasure Hunt

solve a system of ode by bvp4c

3 个评论 显示 1更早的评论隐藏 1更早的评论

采纳的回答

3 个评论 显示 1更早的评论隐藏 1更早的评论

更多回答（1 个）

9 个评论 显示 7更早的评论隐藏 7更早的评论

另请参阅

类别

标签

Community Treasure Hunt

WeChat

3 个评论
显示 1更早的评论隐藏 1更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论

9 个评论
显示 7更早的评论隐藏 7更早的评论