下面是我整段的matlab代码
clear all;
%%程序初始化
drou=0.00001;%密度变化步长
rou1(1)=0.85;%已知初始密度(100-200)
rou2(1)=0.85;%(0- 100)
p1(1)=100;%已知初始压强(100-200)
p2(1)=100;%(0-100)
n=1;%计数器
%%算法主体
while (p1( end)>=0)
n=n+1;
rou1(n)=rou1(n-1) -drou;
i=floor((p1(n-1))/50);
if i==0
E=1538.4;
elseif i==1
E=fun1(p1(n-1));
else
E=fun2(p1(n-1));
end
dp=E/ ((rou1(n-1)+rou1(n))/2)*drou;
p1(n)=p1(n-1)-dp;
end
n=1;
while (p2(end) <=500)
n=n+1;
rou2(n)=rou2(n-1)+drou;
i=floor((p2(n-1)-100)/50)+3;
if i==5
E=3393.4;
elseif i==3
E=fun3(p2(n-1));
else
E=fun4(p2(n-1));
end
dp=E/((rou2(n-1)+rou2(n))/2)*drou;
p2(n)=p2(n-1)+dp;
end
%%数据拼接
for i=1 :size(p1,2)
Result(1,i)=p1(size(p1,2)+1-i);
Result(2,i)=rou1(size(p1,2)+1-i);
end
p2(:,1)=[];
rou2(:,1)=[];
Result=[Resu1t [p2;rou2]];
if Result(1,1) < 0 %去除首尾无效数据
Result(:,1)=[];
end
if Result(1, end)>200
Result(: , end)=[];
end
clear dp drou E i n p1 p2 rou1 rou2;
%%分段多项拟合的函数
function [e]=fun1(p)%0-49.5
e=0.00004*p^3+0.0106*p^2+4.872*p+1538.4;
end
function [e]=fun2(p)%50-99.5
e=0.00006*p^3+0.0062*p^ 2+5.1165* p+1533.7;
end
function [e]=fun3(p)%100-149.5
e=0.0000005*p^4-0.0001*p^3+0.0308*p^2+3.5249*p+1573.3;
end
function [e]=fun4(p) %150- 200
e=0.000001* p^4-0.0007*p^3+0.159*p^2-9.9328*p+2104.3;
end