2

I have 14 first order differential equations. 14 conditions, 7 are initial ones like x1(0)=0, x2(0)=5... 7 are terminal ones x8(10)=25,x9(10)=0....

I think I should use bvp4c

I found this answer but I'm still have a couple of problems:how to solve a system of Ordinary Differential Equations (ODE's) in Matlab

I create a matlab function to put my system in.

   x'=2x
   y'=3x+5y

coding it like:

 xdot=[2x(1);3x(1)+5x(2)]

like I would do in ode45. Then I should do the same for the boundary conditions. But I have no idea on how to code them. I should build a matrix containing them, but I haven't understood how to build it.

I'm trying to use this reference: http://www.math.tamu.edu/~phoward/m442/matode.pdf pag 12, but he does the y2=y' thing and I'm kind lost in my case. also it doesn't explain well how should I put the 14 conditions I have. 7 on one line and 7 on the other? how do I tell the program which value is referred to each variable?

thanks in advance.

here's the actual system. it's a bit huge, so I fear I need numerical methods.

f1=(delta1*gn-(beta*phi*x(7)*x(1)+(1-u1))/(x(1)+x(2)+x(3)+x(4))-mu*x(1)+psi*x(4));
f2=((beta*phi*x(7)*x(1)*(1-u1))/(x(1)+x(2)+x(3)+x(4))-d*x(2)-mu*x(2));
f3=(d*x(2)-(r+r0*u2)*x(3)-(alfa+mu)*x(3));
f4=((r+r0*u2)*x(3)-(mu+phi)*x(4));
f5=(delta2*hp-(phi*teta*x(3)*x(5)*(1-u1))/(x(1)+x(2)+x(3)+x(4))-gamma*x(5));
f6=((phi*teta*x(3)*x(5)*(1-u1))/(x(1)+x(2)+x(3)+x(4))-gamma*x(6)-k*x(6));
f7=(k*x(6)-gamma*x(7));
f8=x(8)*(mu - (beta*phi*x(1)*x(7) - u1 + 1)/(x(1) + x(2) + x(3) + x(4))^2 + (beta*phi*x(7))/(x(1) + x(2) + x(3) + x(4))) + x(9)*((beta*phi*x(7)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4)) - (beta*phi*x(1)*x(7)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2) + (teta*x(12)*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2 - (teta*x(13)*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2;
f9=x(9)*(d + mu - (beta*phi*x(1)*x(7)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2) - d*x(10) - A1 - (x(8)*(beta*phi*x(1)*x(7) - u1 + 1))/(x(1) + x(2) + x(3) + x(4))^2 + (teta*x(12)*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2 - (teta*x(13)*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2;
f10= x(10)*(alfa + mu + r + r0*u2) - A2 - x(11)*(r + r0*u2) - x(12)*((teta*phi*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4)) - (teta*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2) + x(13)*((teta*phi*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4)) - (teta*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2) - (x(8)*(beta*phi*x(1)*x(7) - u1 + 1))/(x(1) + x(2) + x(3) + x(4))^2 - (beta*x(9)*phi*x(1)*x(7)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2;
f11=x(11)*(mu + phi) - x(8)*(psi + (beta*phi*x(1)*x(7) - u1 + 1)/(x(1) + x(2) + x(3) + x(4))^2) - (beta*x(9)*phi*x(1)*x(7)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2 + (teta*x(12)*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2 - (teta*x(13)*phi*x(3)*x(5)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))^2;
f12=x(12)*(gamma - (teta*phi*x(3)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4))) + (teta*x(13)*phi*x(3)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4));
f13=x(13)*(gamma + k) - k*x(14);
f14=gamma*x(14) + (beta*x(8)*phi*x(1))/(x(1) + x(2) + x(3) + x(4)) + (beta*x(9)*phi*x(1)*(u1 - 1))/(x(1) + x(2) + x(3) + x(4));

extra:

u1=max(a1,min(b1,1/(2*B1)*(beta*phi/(x(1)+x(2)+x(3)+x(4))*x(7)*x(1)*(x(9)-x(8))+phi*teta/(x(1)+x(2)+x(3)+x(4))*x(3)*x(5)*(x(13)-x(12)))));
u2=max(a2,min(b2,1/(2*B2)*(r0*x(3)*x(10)-r0*x(3)*x(11))));
Community
  • 1
  • 1
user1834153
  • 330
  • 5
  • 20
  • if your system is just first order, and you have the initial conditions and want numerical solution, you can use ode45. If you have boundary conditions and want numerical solution the use bvp4c, since ode45 is only for initial conditions. You can have a look at http://www.mathworks.com/help/matlab/ref/bvp4c.html for examples. dsolve, since it is symbolic, does not care if the conditions are initial or otherwise. For example, in the solution below, change y(0)==1 to y(1)==1 and it will solve it. – Nasser Mar 20 '15 at 08:40
  • I'm trying dsolve, but I get an error when I add the line regarding u1 and u2. it says they need to be converible to float numbers. it's a bit strange, shound't it take symbolical expressions too? – user1834153 Mar 20 '15 at 08:55
  • I agree that http://se.mathworks.com/help/matlab/ref/bvp4c.html is probably your best best I would be shocked if dsolve can analytically solve a 14 dimensional system, with such nasty non-linearity as max(min()) nested in there – alexandre iolov Mar 20 '15 at 10:09
  • @alexandreiolov thanks, I'm trying bvp4c atm, I have a few problems input the boundary conditions , but this video is helping: https://www.youtube.com/watch?v=iEep1-WnjlM – user1834153 Mar 20 '15 at 10:11

1 Answers1

2

To solve BVP ode's using syms, the ode is y''+3 y' + 3 y = 0, it is first converted to 2 first order (state space formulation) and then solved

clear all; close all
syms x(t) y(t)
Dx  = diff(x);
Dy  = diff(y);
eq1 = Dx == y;
eq2 = Dy == -3*x-5*y;
[x,y] = dsolve(eq1,eq2, x(0) == 0, y(1) ==1)
figure;
ezplot(x,[0,6])

Mathematica graphics

To solve same BVP using bvp4c

clear all
t0 = 0; %initial time
tf = 6; %final time
odefun=@(t,y) [y(2); -3*y(1)-5*y(2)];
bcfun=@(yleft,yright) [yleft(1);yright(1)-1];  
solinit = bvpinit(linspace(t0,tf),[0 1]);

sol = bvp4c(odefun,bcfun,solinit);

figure;
plot(sol.x(1,:),-sol.y(1,:),'r')
title('solution');
xlabel('time');
ylabel('y(t)');
grid;

Mathematica graphics

ps. the numerical solution y-axis value ticks seems not to match the syms. But It looks this is just scaling of value thing. No time to look into it. May be someone can spot something and I'll update.

Nasser
  • 12,849
  • 6
  • 52
  • 104
  • my 14 equations are very huge, I fear I need numerical methods. – user1834153 Mar 20 '15 at 08:23
  • I am using the example you posted. I have no idea what you have. But have you tried them on dsolve first? You said they are first order ode's. I am sure dsolve can solve first order ode's. – Nasser Mar 20 '15 at 08:24
  • if you can have a look at the size of the system and give me an opinion about which kind of solution I shall look for I'd be grateful. – user1834153 Mar 20 '15 at 08:34
  • I tried removing the bound on u1 and u2 but matlab doesn't find an explicit solution. I'll go with bvp4c, could you rewrite your first example with bvp4c syntax and one condition on time=1? thanks in advice. – user1834153 Mar 20 '15 at 09:18
  • it does seem to work on my problem now, just one more question if you don't mind. in plot(sol.x(1,:),-sol.y(1,:),'r') why is there the ''-''? – user1834153 Mar 20 '15 at 17:58
  • @user3149593 The plot I got from bvp4c was, for some reason, flipped relative to what it should be. I do not know why, since I do not use bvp4c much. Might be a configuration issue or option of some sort. Putting the "-" on the y-axis corrected this to make it match the analytical solution. You might want to look at more examples on bvp4c on the net to see if some other setup is needed. – Nasser Mar 20 '15 at 20:12
  • In the numerical solution, the second BC has the wrong index, should be `yright(2)-1` to give the derivative a prescribed value. – Lutz Lehmann Nov 09 '18 at 10:20