1、% 一維對流方程迎風格式、Lax格式、FTCS格式差分法計算% 潭花林 清華大學航天航空學院% FTCS格式對于一維對流方程不穩(wěn)定，最好不用clcclear all % 1.參數(shù)定義dx=1;x1=-18;x2=18;x=x1:dx:x2;L1=length(x);% dt=0.5*dx; % 收斂dt=2*dx; % 不收斂t1=0;t2=t1+80*dt;t=t1:dt:t2;L2=length(t);alpha=1;lambda=alpha*dt/dx;geshi=1; % 迎風格式% geshi=2; % Lax格式% geshi=3; % FTCS格式 % 2.顯式求解zeta=ze

2、ros(L1,L2); for kk=1:3 geshi=kk; for ii=1:L1 if x(ii)0 zeta(ii,1)=1; else if x(ii)=0 zeta(ii,1)=1/2; else if x(ii)0 zeta(ii,1)=0; end end endend if geshi=1 for ii=2:L1 for jj=1:(L2-1) zeta(ii,jj+1)=zeta(ii,jj)-lambda*(zeta(ii,jj)-zeta(ii-1,jj); end zeta(1,jj+1)=zeta(2,jj+1); end zeta1=zeta; else if

3、geshi=2 for ii=2:(L1-1) for jj=1:(L2-1) zeta(ii,jj+1)=(zeta(ii+1,jj)+zeta(ii-1,jj)/2-. lambda/2*(zeta(ii+1,jj)-zeta(ii-1,jj); end zeta(1,jj+1)=zeta(2,jj+1); zeta(L1,jj+1)=zeta(L1,jj)-lambda*(zeta(L1,jj)-zeta(L1-1,jj); end zeta2=zeta;else if geshi=3 for ii=2:(L1-1) for jj=1:(L2-1) zeta(ii,jj+1)=zeta(

4、ii,jj)-lambda/2*(zeta(ii+1,jj)-zeta(ii-1,jj); end zeta(1,jj+1)=zeta(2,jj+1); zeta(L1,jj+1)=zeta(L1,jj)-lambda*(zeta(L1,jj)-zeta(L1-1,jj); end zeta3=zeta; end endendend % 3.繪圖 % 3.1 t=0figure(1)n=1;plot(x,zeta1(1:L1,n),-k,x,zeta2(1:L1,n),-.k,x,zeta3(1:L1,n),-k)% 作圖% axis equal % 是否要求x、y坐標間距相等% grid o

5、n % 是否要求畫網(wǎng)格xlabel(x/m),ylabel(t/s) % x,y軸表示的變量含義%text(1,2,f(x) % 圖中文字標識legend(迎風格式,Lax格式,FTCS格式) % 不同曲線的線型區(qū)分title(t=0時刻的計算結(jié)果) % 標題 axis(-18,18,-0.2,1.2) % 3.2 t=10figure(2)n=(10-t(1)/dt;plot(x,zeta1(1:L1,n),-k,x,zeta2(1:L1,n),-.k,x,zeta3(1:L1,n),-k)% 作圖% axis equal % 是否要求x、y坐標間距相等% grid on % 是否要求畫網(wǎng)格

6、xlabel(x/m),ylabel(t/s) % x,y軸表示的變量含義%text(1,2,f(x) % 圖中文字標識legend(迎風格式,Lax格式,FTCS格式) % 不同曲線的線型區(qū)分title(t=10s時刻的計算結(jié)果) % 標題 % 3.3 t=20figure(3)n=(20-t(1)/dt;plot(x,zeta1(1:L1,n),-k,x,zeta2(1:L1,n),-.k,x,zeta3(1:L1,n),-k)% 作圖% axis equal % 是否要求x、y坐標間距相等% grid on % 是否要求畫網(wǎng)格xlabel(x/m),ylabel(t/s) % x,y軸表示的變量含義%text(1,2,f(x) % 圖中文字標識legend(迎風格式,Lax格式,FTCS格式) % 不同曲線的線型區(qū)分title(t=20s時刻的計算結(jié)果) % 標題 % 3.4 t=40figure(4)n=(40-t(1)/dt;plot(x,zeta1(1:L1,n),-k,x,zeta2(1:L1,n),-.k,x,zeta3(1:L1,n),-k)% 作圖% axis equal % 是否要求x、y坐標間距相等% grid on % 是否要求畫網(wǎng)格xlabel(x/m),yl