如何在 matlab 中解决具有某些终端条件的 ODE 系统？

Question

我有 14 个一阶微分方程。 14 个条件，7 个是初始条件，例如 x1(0)=0、x2(0)=5... 7 是终端的x8(10)=25,x9(10)=0....

我想我应该使用bvp4c

我找到了这个答案，但我还有几个问题：how to solve a system of Ordinary Differential Equations (ODE's) in Matlab

我创建了一个 matlab 函数来放入我的系统。

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

编码如下：

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

就像我在ode45 中所做的那样。 然后我应该对边界条件做同样的事情。但我不知道如何编码它们。 我应该构建一个包含它们的矩阵，但我不知道如何构建它。

我正在尝试使用此参考：http://www.math.tamu.edu/~phoward/m442/matode.pdf 第 12 页，但他做了y2=y' 的事情，我有点迷失在我的案例中。它也没有很好地解释我应该如何放置我拥有的 14 个条件。一条7条，另一条7条？如何告诉程序每个变量引用了哪个值？

提前致谢。

这是实际系统。它有点大，所以我担心我需要数值方法。

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

额外的：

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

Answer 1

使用syms求解BVP ode，ode为y''+3 y' + 3 y = 0，先转换为2一阶（状态空间公式）再求解

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

使用 bvp4c 解决相同的 BVP

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;

ps。数值解 y 轴值刻度似乎与符号不匹配。但看起来这只是价值的缩放。没时间研究它。可能有人可以发现一些东西，我会更新。