clear;
%% 程序初始化
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=[Result [p2;rou2]];
if Result(1,1) 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
显示fun2未命名
function [e]=fun2(p) %50-99.5
这个函数的定义同一个m文件，放在调用前面
或者检查下是不是在同一个目录里
p1(n-1)这个是不是在取值范围内。
clear;
%% 程序初始化
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=[Result [p2;rou2]];
if Result(1,1)
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
显示fun2未命名
如下：
function [e]=fun2(p) %50-99.5
这个函数的定义同一个m文件，放在调用前面
或者检查下是不是在同一个目录里
p1(n-1)这个是不是在取值范围内。