陶文铨--数值传热学-第二版-第五章--5-2

精确解：p=[1,5,10];x=0:1/19:1;fori=1:1:3forj=1:1:20y(i,j)=(exp(p(1,i)*19*x(1,j))-1)/(exp(p(1,i)*19)-1);endplot(x,y(i,:));holdon;end由题对中心差分、一阶迎风、混合格式进行模块编程：他们之间可以通用，只需更改ae关于p的函数即可：程序如下:(1)中心差分p=[1,5,10];fori=1:1:3ae=1-0.5*p(1,i);00.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.70.80.91x/L(Φ-ΦL)/(Φ0-ΦL)精确解图像PΔ=1PΔ=5ΡΔ=10aw=p(1,i)+ae;ap=ae+aw;fori=1:1:18forj=1:1:20a(i,j)=0;endendfori=1:1:18j=i;a(i,j)=aw;a(i,j+1)=-ap;a(i,j+2)=ae;endfori=1:1:17n=i+1;form=i:-1:1b(1,1)=a(m,n);a(m,n)=-a(i+1,n)/a(i+1,n)*b(1,1)+a(m,n);a(m,n+1)=-a(i+1,n+1)/a(i+1,n)*b(1,1)+a(m,n+1);a(m,n+2)=-a(i+1,n+2)/a(i+1,n)*b(1,1)+a(m,n+2);endendF(1)=0;F(20)=1;F(19)=(-a(1,20)*F(20)-a(1,1)*F(1))/a(1,19);fori=2:1:18F(i)=(-a(i,20)*F(20)-a(i,19)*F(19))/a(i,i);endx=0:1/19:1;y(1,:)=F;plot(x,y);holdonend（2）一阶迎风更改ae表达式：ifp(1,i)0ae=1+p(1,i);elseae=1;end（3）混合格式更改ae表达式：ifp(1,i)2ae=0;elseif(p(1,i)=2)&&(p(1,i)0)ae=1-p(1,i)/2;elseae=-p(1,i);end00.10.20.30.40.50.60.70.80.91-0.8-0.6-0.4-0.200.20.40.60.81x/L(Φ-ΦL)/(Φ0-ΦL)中心差分图PΔ=1PΔ=5PΔ=1000.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.70.80.91x/L(Φ-ΦL)/(Φ0-ΦL)一阶迎风图PΔ=1PΔ=5PΔ=1000.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.70.80.91x/L（Φ-ΦL）/(Φ0-ΦL)混合格式图PΔ=1PΔ=5PΔ=10（4）quick格式*12111)35(8122121iiiiiiiPPPPP2i*1111)336(8122121iiiiiiPPPPP2i用matlab编quick格式：程序如下：p1=[1,5,10];fori=1:1:3p=p1(1,i);F(1,1)=0;F(1,20)=1;fori=2:1:19F(1,i)=0;endF1(1,1)=0;F1(1,20)=1;fori=2:1:19F1(1,i)=0;endwhile(F1(1,2)==0)||max0.00000001fori=1:1:20F(1,i)=F1(1,i);endA=1;B=1/(2+p);C=(1+p)/(2+p);D(1,2)=p*((6*F(1,2)-3*F(1,1)-3*F(1,3))/8)/(2+p);fori=3:1:19D(1,i)=p*((5*F(1,i)-F(1,i-1)-F(1,i-2)-3*F(1,i+1))/8)/(2+p);endP(1,2)=B;Q(1,2)=C*F(1,1)+D(1,2);fori=3:1:19P(1,i)=B/(A-C*P(1,i-1));Q(1,i)=(D(1,i)+C*Q(1,i-1))/(A-C*P(1,i-1));endfori=19:-1:2F1(1,i)=P(1,i)*F1(1,i+1)+Q(1,i);endmax=abs(F1(1,2)-F(1,2));fori=3:1:19ifabs(F1(1,i)-F(1,i))maxmax=abs(F1(1,i)-F(1,i));endendendx=0:1/19:1;plot(x,F1(1,:));holdon;end00.10.20.30.40.50.60.70.80.91-0.4-0.200.20.40.60.81x/L(Φ-ΦL)/(Φ0-ΦL)quick格式图PΔ=1PΔ=5PΔ=10