目录
一、第1问
本问主要用到了：偏微分方程转差分方程，再通过矩阵转化利用追赶法求解偏微分方程。
1.1 第一问追赶法求解差分方程与绘图MATLAB程序
%第一问追赶法求解差分方程与绘图MATLAB程序 clear; close all; r=xlsread('2018A.xlsx','附件1','M3:M6'); k=xlsread('2018A.xlsx','附件1','D3:D6'); m=xlsread('2018A.xlsx','附件1','K3:K6'); n=xlsread('2018A.xlsx','附件1','N3'); T0=xlsread('2018A.xlsx','附件2','B5403'); t0=xlsread('2018A.xlsx','附件2','A3:A5403'); a=zeros(sum(m),sum(m)); b=zeros(sum(m),1); u=zeros(sum(m)+1,n+1); he=linspace(8,9,100); for z=1:100 h=0.0001; u(:,1)=37; u(1,:)=75; %a矩阵 for i=1:(sum(m)) if i==1 a(i,i)=1+2*r(1); a(i,i+1)=-r(1); end if (1<i)&&(i<m(1)) a(i,i)=1+2*r(1); a(i,i-1)=-r(1); a(i,i+1)=-r(1); end if (i==m(1)) a(i,i)=(k(1)+k(2))/h; a(i,i-1)=-k(1)/h; a(i,i+1)=-k(2)/h; end if (m(1)<i)&&(i<(m(1)+m(2))) a(i,i)=1+2*r(2); a(i,i-1)=-r(2); a(i,i+1)=-r(2); end if i==(m(1)+m(2)) a(i,i)=(k(2)+k(3))/h; a(i,i-1)=-k(2)/h; a(i,i+1)=-k(3)/h; end if ((m(1)+m(2))<i)&&(i<(m(1)+m(2)+m(3))) a(i,i)=1+2*r(3); a(i,i-1)=-r(3); a(i,i+1)=-r(3); end if i==(m(1)+m(2)+m(3)) a(i,i)=(k(3)+k(4))/h; a(i,i-1)=-k(3)/h; a(i,i+1)=-k(4)/h; end if ((m(1)+m(2)+m(3))<i)&&(i<(m(1)+m(2)+m(3)+m(4))) a(i,i)=1+2*r(4); a(i,i-1)=-r(4); a(i,i+1)=-r(4); end if i==(m(1)+m(2)+m(3)+m(4)) a(i,i)=k(4)/h+he(z); a(i,i-1)=-k(4)/h; end end %b矩阵 for s=2:n+1 for t=1:(sum(m)) if t==1 b(t,1)=u(2,s-1)+r(1)*u(1,s); end if t>1 b(t,1)=u(1+t,s-1); end b(m(1),1)=0; b((m(1)+m(2)),1)=0; b((m(1)+m(2)+m(3)),1)=0; b((m(1)+m(2)+m(3)+m(4)),1)=37*he(z); end %追赶法求解 bb=diag(a)'; aa=[0,diag(a,-1)']; c=diag(a,1)'; N=length(bb); L=zeros(N); uu0=0;y0=0;aa(1)=0; L(1)=bb(1)-aa(1)*uu0; y(1)=(b(1)-y0*aa(1))/L(1); uu(1)=c(1)/L(1); for i=2:(N-1) L(i)=bb(i)-aa(i)*uu(i-1); y(i)=(b(i)-y(i-1)*aa(i))/L(i); uu(i)=c(i)/L(i); end L(N)=bb(N)-aa(N)*uu(N-1); y(N)=(b(N)-y(N-1)*aa(N))/L(N); x(N)=y(N); for i=(N-1):-1:1 x(i)=y(i)-uu(i)*x(i+1); end u(2:sum(m)+1,s)=x'; end if u(153,5401)<=T0 hee=he(z); break end end x=0.1:0.1:15.3; t=1:1:5401; surf(t,x,u) shading interp xlabel('时间(s)') ylabel('厚度(mm)') zlabel('温度(℃)') figure %绘图 T1=u(m(1)+1,:); T2=u(m(2)+m(1)+1,:); T3=u(m(3)+m(2)+m(1)+1,:); T4=u(m(4)+m(3)+m(2)+m(1)+1,:); plot(t0,T1,'-','linewidth',1) axis([0 6000 36 80]); hold on plot(t0,T2,'--','linewidth',1) hold on plot(t0,T3,':.','linewidth',0.5) hold on plot(t0,T4,':','linewidth',1.5) legend('第Ⅰ交界处','第Ⅱ交界处','第Ⅲ交界处','第Ⅳ交界处'); xlabel('时间(s)'); ylabel('温度(℃)'); u=u'; xlswrite('problem1.xlsx',u) %% 第一问 he=8与he=9绘图MATLAB程序 t=xlsread('2018A.xlsx','附件2','A3:A5403'); T=xlsread('2018A.xlsx','附件2','B3:B5403'); TA=xlsread('2018A.xlsx','附件2','C3:C5403'); TB=xlsread('2018A.xlsx','附件2','D3:D5403'); plot(t,T,'-','linewidth',1) axis([0 6000 36 52]); hold on plot(t,TA,':','linewidth',1.5) hold on plot(t,TB,'--','linewidth',1) legend('原始数据','he=8','he=9'); xlabel('时间(s)'); ylabel('温度(℃)'); %p=polyfit(t,T,8); %T2=polyval(p,t); %p=vpa(poly2sym(p),2) %plot(t,T2) %% 第一问 he=8.6162绘图MATLAB程序 t=xlsread('2018A.xlsx','附件2','A3:A5403'); T=xlsread('2018A.xlsx','附件2','B3:B5403'); T0=xlsread('2018A.xlsx','附件2','E3:E5403'); plot(t,T,'-','linewidth',1) axis([0 6000 36 52]); hold on plot(t,T0,'--','linewidth',1) legend('原始数据','he=8.6162'); xlabel('时间(s)'); ylabel('温度(℃)');
二、第2问
本问主要用到了：利用枚举法代入第1问求解出的模型寻找第Ⅱ层厚度最优解。
2.1 第2问枚举法求解最优解与绘图MATLAB程序
%第二问枚举法求解最优解与绘图MATLAB程序 clear; close all; r=xlsread('2018A.xlsx','附件1','M3:M6'); k=xlsread('2018A.xlsx','附件1','D3:D6'); n=3600; h=0.0001; m1=6;m3=36;m4=55; m2=linspace(6,250,245); mm2=250; uu2=zeros((m1+250+m3+m4+1),n+1); for z=1:244 m=m1+m2(z)+m3+m4; a=zeros(m,m); b=zeros(m,1); u=zeros(m+1,n+1); he=8.6162; u(:,1)=37; u(1,:)=65; %a矩阵 for i=1:m if i==1 a(i,i)=1+2*r(1); a(i,i+1)=-r(1); end if (1<i)&&(i<m1) a(i,i)=1+2*r(1); a(i,i-1)=-r(1); a(i,i+1)=-r(1); end if (i==m1) a(i,i)=(k(1)+k(2))/h; a(i,i-1)=-k(1)/h; a(i,i+1)=-k(2)/h; end if (m1<i)&&(i<(m1+m2(z))) a(i,i)=1+2*r(2); a(i,i-1)=-r(2); a(i,i+1)=-r(2); end if i==(m1+m2(z)) a(i,i)=(k(2)+k(3))/h; a(i,i-1)=-k(2)/h; a(i,i+1)=-k(3)/h; end if ((m1+m2(z))<i)&&(i<(m1+m2(z)+m3)) a(i,i)=1+2*r(3); a(i,i-1)=-r(3); a(i,i+1)=-r(3); end if i==(m1+m2(z)+m3) a(i,i)=(k(3)+k(4))/h; a(i,i-1)=-k(3)/h; a(i,i+1)=-k(4)/h; end if ((m1+m2(z)+m3)<i)&&(i<(m1+m2(z)+m3+m4)) a(i,i)=1+2*r(4); a(i,i-1)=-r(4); a(i,i+1)=-r(4); end if i==(m1+m2(z)+m3+m4) a(i,i)=k(4)/h+he; a(i,i-1)=-k(4)/h; end end %b矩阵 for s=2:n+1 for t=1:m if t==1 b(t,1)=u(2,s-1)+r(1)*u(1,s); end if t>1 b(t,1)=u(1+t,s-1); end b(m1,1)=0; b((m1+m2(z)),1)=0; b((m1+m2(z)+m3),1)=0; b((m1+m2(z)+m3+m4),1)=37*he; end %追赶法求解 bb=diag(a)'; aa=[0,diag(a,-1)']; c=diag(a,1)'; N=length(bb); L=zeros(N); uu0=0;y0=0;aa(1)=0; L(1)=bb(1)-aa(1)*uu0; y(1)=(b(1)-y0*aa(1))/L(1); uu(1)=c(1)/L(1); for i=2:(N-1) L(i)=bb(i)-aa(i)*uu(i-1); y(i)=(b(i)-y(i-1)*aa(i))/L(i); uu(i)=c(i)/L(i); end L(N)=bb(N)-aa(N)*uu(N-1); y(N)=(b(N)-y(N-1)*aa(N))/L(N); x(N)=y(N); for i=(N-1):-1:1 x(i)=y(i)-uu(i)*x(i+1); end u(2:m+1,s)=x'; end if (u(m+1,3601))<=47&&(u(m+1,3301)<44) if mm2>m2(z) mm2=m2(z); uu2=u; end end end %% 第二问 第Ⅱ层厚度=19.2mm绘图MATLAB程序 t=xlsread('2018A.xlsx','附件2','A3:A3603'); T2=xlsread('2018A.xlsx','附件2','J3:J3603'); plot(t,T2,'-','linewidth',1) grid on legend('第Ⅱ层厚度=19.2mm'); xlabel('时间(s)'); ylabel('温度(℃)'); %% 第二问 灵敏度分析MATLAB程序 T=xlsread('2018A.xlsx','附件2','G1:L1'); d=xlsread('2018A.xlsx','附件2','G2:L2'); plot(T,d,'linewidth',1) xlabel('温度(℃)'); ylabel('第Ⅱ层厚度(mm)');
三、第3问
本问主要用到了：继续利用枚举法，通过for循环嵌套，代入第1问求解出的模型寻找第Ⅱ层与第Ⅳ层厚度的最优解。并通过比较“以厚度为衡量标准”与“以重量为衡量标准”，最终确定“以重量为衡量标准”，可以在同等厚度的情况下，使得服装的重量更轻。
3.1 第三问枚举法求解最优解与绘图MATLAB程序
%第三问枚举法求解最优解与绘图MATLAB程序（以厚度为衡量标准） clear; close all; r=xlsread('2018A.xlsx','附件1','M3:M6'); k=xlsread('2018A.xlsx','附件1','D3:D6'); n=1800; h=0.0001; m1=6;m3=36; m2=linspace(6,250,245); m4=linspace(6,64,59); mm2=250;mm4=64; uu2=zeros((m1+250+m3+64+1),n+1); for v=1:58 clear x; for z=1:244 m=m1+m2(z)+m3+m4(v); a=zeros(m,m); b=zeros(m,1); u=zeros(m+1,n+1); he=8.6162; u(:,1)=37; u(1,:)=80; %a矩阵 for i=1:m if i==1 a(i,i)=1+2*r(1); a(i,i+1)=-r(1); end if (1<i)&&(i<m1) a(i,i)=1+2*r(1); a(i,i-1)=-r(1); a(i,i+1)=-r(1); end if (i==m1) a(i,i)=(k(1)+k(2))/h; a(i,i-1)=-k(1)/h; a(i,i+1)=-k(2)/h; end if (m1<i)&&(i<(m1+m2(z))) a(i,i)=1+2*r(2); a(i,i-1)=-r(2); a(i,i+1)=-r(2); end if i==(m1+m2(z)) a(i,i)=(k(2)+k(3))/h; a(i,i-1)=-k(2)/h; a(i,i+1)=-k(3)/h; end if ((m1+m2(z))<i)&&(i<(m1+m2(z)+m3)) a(i,i)=1+2*r(3); a(i,i-1)=-r(3); a(i,i+1)=-r(3); end if i==(m1+m2(z)+m3) a(i,i)=(k(3)+k(4))/h; a(i,i-1)=-k(3)/h; a(i,i+1)=-k(4)/h; end if ((m1+m2(z)+m3)<i)&&(i<(m1+m2(z)+m3+m4(v))) a(i,i)=1+2*r(4); a(i,i-1)=-r(4); a(i,i+1)=-r(4); end if i==(m1+m2(z)+m3+m4(v)) a(i,i)=k(4)/h+he; a(i,i-1)=-k(4)/h; end end %b矩阵 for s=2:n+1 for t=1:m if t==1 b(t,1)=u(2,s-1)+r(1)*u(1,s); end if t>1 b(t,1)=u(1+t,s-1); end b(m1,1)=0; b((m1+m2(z)),1)=0; b((m1+m2(z)+m3),1)=0; b((m1+m2(z)+m3+m4(v)),1)=37*he; end %追赶法求解 bb=diag(a)'; aa=[0,diag(a,-1)']; c=diag(a,1)'; N=length(bb); L=zeros(N); uu0=0;y0=0;aa(1)=0; L(1)=bb(1)-aa(1)*uu0; y(1)=(b(1)-y0*aa(1))/L(1); uu(1)=c(1)/L(1); for i=2:(N-1) L(i)=bb(i)-aa(i)*uu(i-1); y(i)=(b(i)-y(i-1)*aa(i))/L(i); uu(i)=c(i)/L(i); end L(N)=bb(N)-aa(N)*uu(N-1); y(N)=(b(N)-y(N-1)*aa(N))/L(N); x(N)=y(N); for i=(N-1):-1:1 x(i)=y(i)-uu(i)*x(i+1); end u(2:m+1,s)=x'; end if (u(m+1,1801))<=47&&(u(m+1,1501)<44) if (mm2+mm4)>(m2(z)+m4(v)) mm2=m2(z); mm4=m4(v); uu2=u; end end end end %第三问枚举法求解最优解与绘图MATLAB程序（以重量为衡量标准） clear; close all; r=xlsread('2018A.xlsx','附件1','M3:M6'); k=xlsread('2018A.xlsx','附件1','D3:D6'); n=1800; h=0.0001; m1=6;m3=36; m2=linspace(6,250,245); m4=linspace(6,64,59); mm2=250;mm4=64; uu2=zeros((m1+250+m3+64+1),n+1); for v=1:58 clear x; for z=1:244 m=m1+m2(z)+m3+m4(v); a=zeros(m,m); b=zeros(m,1); u=zeros(m+1,n+1); he=8.6162; u(:,1)=37; u(1,:)=80; %a矩阵 for i=1:m if i==1 a(i,i)=1+2*r(1); a(i,i+1)=-r(1); end if (1<i)&&(i<m1) a(i,i)=1+2*r(1); a(i,i-1)=-r(1); a(i,i+1)=-r(1); end if (i==m1) a(i,i)=(k(1)+k(2))/h; a(i,i-1)=-k(1)/h; a(i,i+1)=-k(2)/h; end if (m1<i)&&(i<(m1+m2(z))) a(i,i)=1+2*r(2); a(i,i-1)=-r(2); a(i,i+1)=-r(2); end if i==(m1+m2(z)) a(i,i)=(k(2)+k(3))/h; a(i,i-1)=-k(2)/h; a(i,i+1)=-k(3)/h; end if ((m1+m2(z))<i)&&(i<(m1+m2(z)+m3)) a(i,i)=1+2*r(3); a(i,i-1)=-r(3); a(i,i+1)=-r(3); end if i==(m1+m2(z)+m3) a(i,i)=(k(3)+k(4))/h; a(i,i-1)=-k(3)/h; a(i,i+1)=-k(4)/h; end if ((m1+m2(z)+m3)<i)&&(i<(m1+m2(z)+m3+m4(v))) a(i,i)=1+2*r(4); a(i,i-1)=-r(4); a(i,i+1)=-r(4); end if i==(m1+m2(z)+m3+m4(v)) a(i,i)=k(4)/h+he; a(i,i-1)=-k(4)/h; end end %b矩阵 for s=2:n+1 for t=1:m if t==1 b(t,1)=u(2,s-1)+r(1)*u(1,s); end if t>1 b(t,1)=u(1+t,s-1); end b(m1,1)=0; b((m1+m2(z)),1)=0; b((m1+m2(z)+m3),1)=0; b((m1+m2(z)+m3+m4(v)),1)=37*he; end %追赶法求解 bb=diag(a)'; aa=[0,diag(a,-1)']; c=diag(a,1)'; N=length(bb); L=zeros(N); uu0=0;y0=0;aa(1)=0; L(1)=bb(1)-aa(1)*uu0; y(1)=(b(1)-y0*aa(1))/L(1); uu(1)=c(1)/L(1); for i=2:(N-1) L(i)=bb(i)-aa(i)*uu(i-1); y(i)=(b(i)-y(i-1)*aa(i))/L(i); uu(i)=c(i)/L(i); end L(N)=bb(N)-aa(N)*uu(N-1); y(N)=(b(N)-y(N-1)*aa(N))/L(N); x(N)=y(N); for i=(N-1):-1:1 x(i)=y(i)-uu(i)*x(i+1); end u(2:m+1,s)=x'; end if (u(m+1,1801))<=47&&(u(m+1,1501)<44) if (862*mm2+1.18*mm4)>(862*m2(z)+1.18*m4(v)) mm2=m2(z); mm4=m4(v); uu2=u; end end end end %% 第三问 以厚度为衡量标准绘图MATLAB程序 t=xlsread('2018A.xlsx','附件2','A3:A1803'); T3=xlsread('2018A.xlsx','附件2','S3:S1803'); plot(t,T3,'-','linewidth',1) axis([0 2000 36 46]); grid on legend('第Ⅱ层厚度=21.8mm，第Ⅳ层厚度=6.2mm'); xlabel('时间(s)'); ylabel('温度(℃)'); %% 第三问 以重量为衡量标准绘图MATLAB程序 t=xlsread('2018A.xlsx','附件2','A3:A1803'); T3=xlsread('2018A.xlsx','附件2','T3:T1803'); plot(t,T3,'-','linewidth',1) axis([0 2000 36 46]); grid on legend('第Ⅱ层厚度=21.7mm，第Ⅳ层厚度=6.3mm'); xlabel('时间(s)'); ylabel('温度(℃)'); %% 第三问 灵敏度分析MATLAB程序 T=xlsread('2018A.xlsx','附件2','V1:AB1'); d2=xlsread('2018A.xlsx','附件2','V2:AB2'); d4=xlsread('2018A.xlsx','附件2','V3:AB3'); plot(T,d2,'linewidth',1) xlabel('温度(℃)'); ylabel('第Ⅱ层厚度(mm)'); figure plot(T,d4,'linewidth',1) xlabel('温度(℃)'); ylabel('第Ⅳ层厚度(mm)');
355
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
