我正在尝试编写MatLab代码,该代码可以扩展要解决的ODE数量。这是我目前拥有的代码,为简单起见,我从弹簧质量阻尼器系统开始。
clear;clc;close all;
tspan=[0:0.01:1];
x0=[1;0;-1;0;7;0;-7;0;5;0;-5;0;10;0;-10;0];
[t,x] = ode45(@Spring_Mass_Damper,tspan,x0);
figure(1)
plot(t,x(:,1));
hold on
plot(t,3));
hold on
plot(t,5));
hold on
plot(t,7));
hold on
plot(t,9));
hold on
plot(t,11));
hold on
plot(t,13));
hold on
plot(t,15));
grid on
xlabel('Time')
ylabel('Displacement(x)')
title('Displacement Vs Time')
function xp = Spring_Mass_Damper(t,x)
c1=10;
G=9.81;
m1=1;
k1=2000;
xp=[[x(2)];(-(c1/m1).*[x(2)]-(k1/m1).*[x(1)])-(G.*m1);
[x(4)];(-(c1/m1).*[x(4)]-(k1/m1).*[x(3)])-(G.*m1);
[x(6)];(-(c1/m1).*[x(6)]-(k1/m1).*[x(5)])-(G.*m1);
[x(8)];(-(c1/m1).*[x(8)]-(k1/m1).*[x(7)])-(G.*m1);
[x(10)];(-(c1/m1).*[x(10)]-(k1/m1).*[x(9)])-(G.*m1);
[x(12)];(-(c1/m1).*[x(12)]-(k1/m1).*[x(11)])-(G.*m1);
[x(14)];(-(c1/m1).*[x(14)]-(k1/m1).*[x(13)])-(G.*m1);
[x(16)];(-(c1/m1).*[x(16)]-(k1/m1).*[x(15)])-(G.*m1)];
end
我所关心的主要问题是代码的这一领域:
xp=[[x(2)];(-(c1/m1).*[x(2)]-(k1/m1).*[x(1)])-(G.*m1);
[x(4)];(-(c1/m1).*[x(4)]-(k1/m1).*[x(3)])-(G.*m1);
[x(6)];(-(c1/m1).*[x(6)]-(k1/m1).*[x(5)])-(G.*m1);
[x(8)];(-(c1/m1).*[x(8)]-(k1/m1).*[x(7)])-(G.*m1);
[x(10)];(-(c1/m1).*[x(10)]-(k1/m1).*[x(9)])-(G.*m1);
[x(12)];(-(c1/m1).*[x(12)]-(k1/m1).*[x(11)])-(G.*m1);
[x(14)];(-(c1/m1).*[x(14)]-(k1/m1).*[x(13)])-(G.*m1);
[x(16)];(-(c1/m1).*[x(16)]-(k1/m1).*[x(15)])-(G.*m1)];
有没有一种简单的方法可以实现相同数量的系统(ODE),而无需在块中复制和粘贴第一个方程式并手动更改x索引?