function main()
%输入sishuixiangrk4('f1','f2','f3','f4',0,10,12,12.5,1.6,1.5,500)获取结果
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
%x'=f(t,x,y) y'=g(t,x,y)
%N为迭代次数
%h为步长
%ya,xa为初值
in1=[];
in2=[];
in3=[];
in4=[];
in1=3.33.*(0.7+normrnd(0,0.01,1,501)).*(3+normrnd(0,0.05,1,501))/28;
in2=3.33.*(0.6+normrnd(0,0.01,1,501)).*(3+normrnd(0,0.05,1,501))/32;
in3=3.33.*(0.4+normrnd(0,0.01,1,501)).*(3+normrnd(0,0.05,1,501))/28;
in4=3.33.*(0.3+normrnd(0,0.01,1,501)).*(3+normrnd(0,0.05,1,501))/32;
fs1=in1*28;
fs2=in2*32;
fs3=in3*28;
fs4=in4*32;
for i=1:501;
f1=@(t,h1,h3)(-0.11233*sqrt(h1)+0.11233*sqrt(h3)+in1(:,i));
f2=@(t,h2,h4)(-0.07884*sqrt(h2)+0.07884*sqrt(h4)+in2(:,i));
f3=@(t,h1,h3)(-0.11233*sqrt(h3)+in3(:,i));
f4=@(t,h2,h4)(-0.07884*sqrt(h4)+in4(:,i));
h=0.1;
T=zeros(1,501);
h1=zeros(1,501);
h2=zeros(1,501);
h3=zeros(1,501);
h4=zeros(1,501);
T=0:0.1:50;
h1(1)=12.3;
h2(1)=12.65;
h3(1)=1.6;
h4(1)=1.45;
for j=1:500
f11=feval(f1,T(j),h1(j),h3(j));
f31=feval(f3,T(j),h1(j),h3(j));
f12=feval(f1,T(j)+h/2,h1(j)+h/2*f11,h3(j)+f31*h/2);
f32=feval(f3,T(j)+h/2,h1(j)+h/2*f11,h3(j)+h/2*f31);
f13=feval(f1,T(j)+h/2,h1(j)+h/2*f12,h3(j)+h*f32/2);
f33=feval(f3,T(j)+h/2,h1(j)+h/2*f12,h3(j)+h/2*f32);
f14=feval(f1,T(j)+h,h1(j)+h*f13,h3(j)+h*f33);
f34=feval(f3,T(j)+h,h1(j)+h*f13,h3(j)+h*f33);
h1(j+1)=h1(j)+h*(f11+2*f12+2*f13+f14)/6;
h3(j+1)=h3(j)+h*(f31+2*f32+2*f33+f34)/6;
f21=feval(f2,T(j),h2(j),h4(j));
f41=feval(f4,T(j),h2(j),h4(j));
f22=feval(f2,T(j)+h/2,h2(j)+h/2*f21,h4(j)+f41*h/2);
f42=feval(f4,T(j)+h/2,h2(j)+h/2*f21,h4(j)+h/2*f41);
f23=feval(f2,T(j)+h/2,h2(j)+h/2*f22,h4(j)+h*f42/2);
f43=feval(f4,T(j)+h/2,h2(j)+h/2*f22,h4(j)+h/2*f42);
f24=feval(f2,T(j)+h,h2(j)+h*f23,h4(j)+h*f43);
f44=feval(f4,T(j)+h,h2(j)+h*f23,h4(j)+h*f43);
h2(j+1)=h2(j)+h*(f21+2*f22+2*f23+f24)/6;
h4(j+1)=h4(j)+h*(f41+2*f42+2*f43+f44)/6;
R=[(h1+normrnd(0,0.1,1,501))' (h2+normrnd(0,0.1,1,501))' (h3+normrnd(0,0.1,1,501))' (h4+normrnd(0,0.1,1,501))'];%501*4
I=[in1' in2' in3' in4'];%501*4
end
end
h4f=zeros(501,1);
for i=201:501
h4f(i,1)=0.8;
end
h4f
Rf=[R(:,1) R(:,2) R(:,3) R(:,4)+h4f];
Xf=[Rf I];%501*8
Xf
figure(1)
subplot(2,2,1)
plot(R(:,1))
title('水箱1液位')
axis([0,500,11.5,13])
subplot(2,2,2)
plot(R(:,2))
title('水箱2液位')
axis([0,500,12,13])
subplot(2,2,3)
plot(R(:,3))
title('水箱3液位')
axis([0,500,1.2,2])
subplot(2,2,4)
plot(R(:,4))
title('水箱4液位')
axis([0,500,1,2])
figure(2)
subplot(2,2,1)
plot(fs1)
title('水流率f1')
axis([0,500,6,7.5])
subplot(2,2,2)
plot(fs2)
title('水流率f2')
axis([0,500,5.5,6.5])
subplot(2,2,3)
plot(fs3)
title('水流率f3')
axis([0,500,3.5,4.5])
subplot(2,2,4)
plot(fs4)
title('水流率f4')
axis([0,500,2.5,3.5])
save Xf.dat Xf -ascii
- 1
- 2
前往页