24小时热门版块排行榜

返回列表

chufufang

金虫 (小有名气)

应助: 2 (幼儿园)
金币: 462.8
红花: 2
帖子: 244
在线: 119.3小时
虫号: 2384076
注册: 2013-03-28
性别: MM
专业: 控制理论与方法

[交流] 离散微分方程代码已有1人参与

clc,clear
n=100;
x1(1)=1;
x2(1)=3;
x3(1)=5;
x4(1)=7;
x5(1)=9;
x6(1)=8;
x7(1)=6;
x8(1)=4;

x9(1)=0.2;
x10(1)=0.4;
x11(1)=0.6;
x12(1)=0.8;
x13(1)=0.9;
x14(1)=0.7;
x15(1)=0.5;
x16(1)=0.3;

k(1)=0;

for i=2:n;
T=0.3;
c=1;
d=2;
c1=1;
x1(i)=x1(i-1)+T*x9(i-1)-c*T*T*(((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2-d^2)/((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2)^2)*x1(i-1)+c*T*T*(((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2-d^2)/((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2)^2)*x2(i-1)-1/2*T*T*c1*x9(i-1)+1/2*T*T*c1*x10(i-1);
x2(i)=x2(i-1)+T*x10(i-1)-c*T*T*(((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2-d^2)/((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2)^2)*x2(i-1)+c*T*T*(((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2-d^2)/((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2)^2)*x3(i-1)-1/2*T*T*c1*x10(i-1)+1/2*T*T*c1*x11(i-1);
x3(i)=x3(i-1)+T*x11(i-1)-c*T*T*(((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2-d^2)/((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2)^2)*x3(i-1)+c*T*T*(((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2-d^2)/((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2)^2)*x4(i-1)-1/2*T*T*c1*x11(i-1)+1/2*T*T*c1*x12(i-1);
x4(i)=x4(i-1)+T*x12(i-1)-c*T*T*(((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2-d^2)/((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2)^2)*x4(i-1)+c*T*T*(((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2-d^2)/((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2)^2)*x1(i-1)-1/2*T*T*c1*x12(i-1)+1/2*T*T*c1*x9(i-1);

x5(i)=x5(i-1)+T*x13(i-1)-c*T*T*(((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2-d^2)/((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2)^2)*x5(i-1)+c*T*T*(((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2-d^2)/((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2)^2)*x6(i-1)-1/2*T*T*c1*x13(i-1)+1/2*T*T*c1*x14(i-1);
x6(i)=x6(i-1)+T*x14(i-1)-c*T*T*(((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2-d^2)/((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2)^2)*x6(i-1)+c*T*T*(((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2-d^2)/((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2)^2)*x7(i-1)-1/2*T*T*c1*x14(i-1)+1/2*T*T*c1*x15(i-1);
x7(i)=x7(i-1)+T*x15(i-1)-c*T*T*(((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2-d^2)/((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2)^2)*x7(i-1)+c*T*T*(((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2-d^2)/((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2)^2)*x8(i-1)-1/2*T*T*c1*x15(i-1)+1/2*T*T*c1*x16(i-1);
x8(i)=x8(i-1)+T*x16(i-1)-c*T*T*(((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2-d^2)/((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2)^2)*x8(i-1)+c*T*T*(((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2-d^2)/((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2)^2)*x5(i-1)-1/2*T*T*c1*x16(i-1)+1/2*T*T*c1*x13(i-1);

x9(i)=x9(i-1)-2*T*c*(((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2-d^2)/((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2)^2)*x1(i-1)+2*c*T*(((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2-d^2)/((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2)^2)*x2(i-1)-T*c1*x9(i-1)+T*c1*x10(i-1);
x10(i)=x10(i-1)-2*T*c*(((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2-d^2)/((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2)^2)*x2(i-1)+2*c*T*(((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2-d^2)/((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2)^2)*x3(i-1)-T*c1*x10(i-1)+T*c1*x11(i-1);
x11(i)=x11(i-1)-2*T*c*(((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2-d^2)/((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2)^2)*x3(i-1)+2*c*T*(((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2-d^2)/((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2)^2)*x4(i-1)-T*c1*x11(i-1)+T*c1*x12(i-1);
x12(i)=x12(i-1)-2*T*c*(((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2-d^2)/((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2)^2)*x4(i-1)+2*c*T*(((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2-d^2)/((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2)^2)*x1(i-1)-T*c1*x12(i-1)+T*c1*x9(i-1);

x13(i)=x13(i-1)-2*T*c*(((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2-d^2)/((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2)^2)*x5(i-1)+2*c*T*(((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2-d^2)/((x2(i-1)-x1(i-1))^2+(x6(i-1)-x5(i-1))^2)^2)*x6(i-1)-T*c1*x13(i-1)+T*c1*x14(i-1);
x14(i)=x14(i-1)-2*T*c*(((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2-d^2)/((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2)^2)*x6(i-1)+2*c*T*(((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2-d^2)/((x3(i-1)-x2(i-1))^2+(x7(i-1)-x6(i-1))^2)^2)*x7(i-1)-T*c1*x14(i-1)+T*c1*x15(i-1);
x15(i)=x15(i-1)-2*T*c*(((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2-d^2)/((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2)^2)*x7(i-1)+2*c*T*(((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2-d^2)/((x4(i-1)-x3(i-1))^2+(x8(i-1)-x7(i-1))^2)^2)*x8(i-1)-T*c1*x15(i-1)+T*c1*x16(i-1);
x16(i)=x16(i-1)-2*T*c*(((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2-d^2)/((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2)^2)*x8(i-1)+2*c*T*(((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2-d^2)/((x1(i-1)-x4(i-1))^2+(x5(i-1)-x8(i-1))^2)^2)*x5(i-1)-T*c1*x16(i-1)+T*c1*x13(i-1);
k(i)=k(i-1)+1;
end

i=1:n;
% plot(k,x9(