如何在 matlab 中解决具有某些终端条件的 ODE 系统?
how can I solve a system of ODEs with some terminal conditions in matlab?
我有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))));
求解BVP ode's using syms,ode为y''+3 y' + 3 y = 0
,先转为2一阶(状态space公式),然后求解
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 轴值刻度似乎与符号不匹配。但看起来这只是价值的缩放。没有时间研究它。可能有人能发现一些东西,我会更新。
我有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))));
求解BVP ode's using syms,ode为y''+3 y' + 3 y = 0
,先转为2一阶(状态space公式),然后求解
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 轴值刻度似乎与符号不匹配。但看起来这只是价值的缩放。没有时间研究它。可能有人能发现一些东西,我会更新。