19年全国数学建模比赛A题代码（简单的迭代思想）

前言：怀念通宵刷题的日子，困得一批

• 第一题

M = 10000; %循环次数
%需要用到的函数
%function E = EMPA(p) 计算弹性模量 p为当前时刻压强
%funtion QA = calcQA(Pnow,rouA,pA) 计算下一时刻A口的流速

%定义好的量,不会参与迭代的量
V = 500 * pi * (10/2) ^ 2;
pA = 160; % MPA 高压油管的压强
tubep0 = 100; % MPA 管子内初始时刻的压强
tao = 1e-3; %迭代的关键点，控制步长
filename = 'myfile.xlsx';
RRange = 'D1';
string0 = 'D';
sheet = 1;
%t0时刻
Et0 = EMPA(tubep0);
rou0 = 0.85;
rouA = rou0 + (pA - tubep0) / Et0 * rou0;
QA0 = calcQA(tubep0,rouA,pA);

tubep = tubep0;
rou = rou0;
num = 0;
while (tubep - tubep0) < 1e-1
    Et = EMPA(tubep);
    QA = calcQA(tubep,rouA,pA);
    deltaROUtao = tao * QA * rouA / V;
    deltaPtao = Et / rou * deltaROUtao;
    rou = rou + deltaROUtao;
    tubep = tubep + deltaPtao;
    %fprintf('the P/MPA of tube is%.6f\n',tubep)
    num = num + 1;
    RRange = strcat(string0,num2str(num));
    xlswrite(filename,tubep,sheet,RRange);
    if (num >= M)
    disp('Warnning:迭代次数太多，可能不收敛');
    return;
    end
end

给出函数的定义

function QA = calcQA(Pnow,rouA,PA)
    QA = 0.85 * (1.4/2)^2 * pi * sqrt(2*(PA-Pnow)/rouA);
end

function E = EMPA(p0)
    E = (0.000000100037752*p0^3-0.000001082481397*p0^2+0.005474444341310*p0+1.531868405848777)*1e3;   
end

• A和B混合注入的情况

filename = 'myfile.xlsx';
RRange = 'A1';
sheet = 3;
string0 = 'A';
%先从起点开始AB混合作用

V = 500 * pi * (10/2)^2;

%定义初始时刻
tubePat0 = 100;
tao = 1e-2;
BOMB = 1e4;
num = 0;
PaA = 160;% MPA 高压管的压强

%t0时刻
Et0 = EMPA(tubePat0);
rouTube = 0.85;
rouA = rouTube + (PaA - tubePat0) / Et0 * rouTube;
rouB = rouTube;

tubeP = tubePat0;
rou = rouTube;
timeB = 0;
while (tubeP - tubePat0) < 5e-1
    Et = EMPA(tubeP);
    QA = calcQA(tubeP,rouA,PaA);
    QB = tixingB(timeB);
    deltaROUtao = tao * (QA * rouA + QB * rouB) / V;
    deltaPtao = Et * deltaROUtao / rou;
    rou = rou + deltaROUtao;
    tubeP = tubeP + deltaPtao;
    num = num + 1;
    tubecancha = (tubeP - tubePat0) * 1e3;
    RRange = strcat(string0,num2str(num));
    xlswrite(filename,tubecancha,sheet,RRange);
    timeB = num * tao;
    if (num >= BOMB)
        disp('Warning:迭代次数太多，可能不收敛');
        return;
    end
end

function QB = tixingB(timeB)
    if timeB < 0.2
        QB = (-1) * 100 * timeB;
    elseif timeB < 2.2
        QB = (-1) * 20;
    elseif timeB < 2.4
        QB = (-1) * (20 - 100 * (timeB - 2.2));
    else
        QB = 0;
    end
end

• 第三题

V = 500 * pi * (10/2)^2;
r1 = 7.239;
r2 = 2.413;
rou  = 0.85 + (0.5 - 100) / (EMPA(100) / 0.85);
Vt0 = (r1-r2)*(5/2)^2*pi + 20;
M  = Vt0 * rou;

ptube = 100;
routube = 0.85;
tao = 0.1;
num = t / tao;
sign = 0;
numtemp = 0;
while( numtemp < num )
	theta = numtemp *tao* w + 3.14;
	rt = calcJiJing(theta);
	Vt = (r1 - rt)*(5/2)^2*pi+20;
	rout = M / Vt;
	Pt = (rout -0.85)/0.85*EMPA(100)+100;
	if(Pt > ptube)
		if((Pt - ptube > 1) && sign == 0)
		%开始放洪
            sign = 1;
            ptube = ptube - 2;
        end
	Et = EMPA(ptube);
	QA = 0.85*(1.4/2)^2*pi*sqrt(2*(Pt-ptube)/rout);
	deltarou = tao*QA*rout/V;
	deltaP = Et/routube*deltarou;
	routube = routube + deltarou;
	ptube = ptube + deltaP;
	M = M - QA *tao;
	end
	numtemp = numtemp + 1;
    RRange = strcat('B',num2str(numtemp));
    xlswrite('try3.xlsx',Pt,4,RRange)
end

function JiJing = calcJiJing(JiJiao)
    JiJing = 0.002396583071159*JiJiao^6-0.045169447134519*JiJiao^5+0.259343667669753*JiJiao^4 -0.287191167670985*JiJiao^3-0.949901781277680*JiJiao^2-0.094085519017525*JiJiao+ 7.247157138979478;
end

• 总结

下列代码可以防止程序陷入死循环导致电脑死机

%把这个放入循环作为循环结束的条件
%if (num >= M)
%    disp('Warnning:迭代次数太多，可能不收敛');
%    return;
%end

• 惯用手法：最小二乘法拟合

filename = ('myfile.xlsx');
x = xlsread(filename,'A2:A402');
y = xlsread(filename,'B2:B402');
%y1 = log(y);
A = polyfit(x,y,3);
X1 = 0:0.01:10;
p1 = polyval(A,X1);
plot(X1,p1,x,y,'o')