NHT.zip

clc clear all; n=0; x=[-10:0.05:10]; for i=-10:0.05:10%初始点坐标的分布 n=n+1; if(i>=-0.5&&i<0) Wave_U(1,n)=sqrt(3)*i+sqrt(3)/2; %坐标点的y方向坐标 elseif (i>=0&&i<=0.5) Wave_U(1,n)=-sqrt(3)*i+sqrt(3)/2;%坐标点的y方向坐标 else Wave_U(1,n)=0; end if(i<=0) Shock_U(1,n)=1.0; else Shock_U(1,n)=0.0; end if(i<=0) rarefaction(1,n)=0.0; else rarefaction(1,n)=1.0; end end diagram=input("select initial data:triangle(1),shock(2),rarefaction(3)"); switch diagram %选择不同的初始点 case(1) equation=input("select equation:wave(1),heat(2),burger(3)"); %选择要计算的方程 switch equation case(1) DT=0.02;%Δt C=2;%c method=input("select method:\nEuler explicit Mehods(1),Upstream Differencing Method(2),Second Order Upwind Method(3),\nLeap Frog(4),Lax-Wendroff(5),MacCormac(6)"); %选择有限差分的方法 switch method case(1)%Euler explicit Mehods for t=2:200 for i=1:n if(i>=2&&i<=n-1) Wave_U(t,i)=Wave_U(t-1,i)-(C*DT/0.05)*(Wave_U(t-1,i+1)-Wave_U(t-1,i));%wave equation_Euler Explicit Method elseif(i==1) Wave_U(t,1)=Wave_U(t-1,1); else Wave_U(t,n)=Wave_U(t-1,n); end x(t,i)=x(t-1,i); end end case(2)%Upstream Differencing Method for t=2:200 for i=1:n if(i>=2&&i<=n-1) Wave_U(t,i)=Wave_U(t-1,i)-(C*DT/0.05)*(Wave_U(t-1,i)-Wave_U(t-1,i-1));%Upstream Differencing Method elseif(i==1) Wave_U(t,1)=Wave_U(t-1,1); else Wave_U(t,n)=Wave_U(t-1,n); end x(t,i)=x(t-1,i); end end case(3)%Second Order Upwind Method for t=2:200 for i=1:n if(i>=2&&i<=n-1) Ave_U(t,i)=Wave_U(t-1,i)-(C*DT/0.05)*(Wave_U(t-1,i)-Wave_U(t-1,i-1));%wave equation_Upstream Differencing Method Wave_U(t,i)=0.5*(Wave_U(t-1,i)+Ave_U(t,i)-(C*DT/0.05)*(Ave_U(t,i)-Ave_U(t,i-1))-(C*DT/0.05)*(Wave_U(t-1,i)-Wave_U(t-1,i-1))); elseif(i==1) Ave_U(t,1)=Wave_U(t-1,1); Wave_U(t,1)=Wave_U(t-1,1); else Ave_U(t,n)=Wave_U(t-1,n); Wave_U(t,n)=Wave_U(t-1,n); end x(t,i)=x(t-1,i); end end case(4)%Leap Frog Method for t=2:200 for i=1:n if(i>=2&&i<=n-1) Wave_U(t,i)=Wave_U(t-1,i)-(C*DT/0.05)*(Wave_U(t-1,i+1)-Wave_U(t-1,i-1)); elseif(i==1) Wave_U(t,1)=Wave_U(t-1,1); else Wave_U(t,n)=Wave_U(t-1,n); end x(t,i)=x(t-1,i); end end case(5)%Lax-Wendroff Method for t=2:200 for i=1:n if(i>=2&&i<=n-1) Wave_U(t,i)=Wave_U(t-1,i)-(0.5*C*DT/0.05)*(Wave_U(t-1,i+1)-Wave_U(t-1,i))+(2*(0.5*C*DT/0.05)^2)*(Wave_U(t-1,i+1)-2*Wave_U(t-1,i)+Wave_U(t-1,i-1)); elseif(i==1) Wave_U(t,1)=Wave_U(t-1,1); else Wave_U(t,n)=Wave_U(t-1,n); end x(t,i)=x(t-1,i); end end case(6)% MacCormac Method for t=2:200 for i=1:n if(i>=2&&i<=n-1) Ave_U(t,i)=Wave_U(t-1,i)-(C*DT/0.05)*(Wave_U(t-1,i+1)-Wave_U(t-1,i)); Wave_U(t,i)=0.5*(Wave_U(t-1,i)+Ave_U(t,i)-(C*DT/0.05)*(Ave_U(t,i)-Ave_U(t,i-1))); elseif(i==1) Ave_U(t,1)=Wave_U(t-1,1)-(C*DT/0.05)*(Wave_U(t-1,2)-Wave_U(t-1,1)); Wave_U(t,1)=Wave_U(t-1,1); else Ave_U(t,n)=Wave_U(t-1,n); Wave_U(t,n)=Wave_U(t-1,n); end x(t,i)=x(t-1,i); end end end %画图 axis([-10 10 -1 5]) plot(x(1,:),Wave_U(1,:),'-r') hold on for i=10:20:200 plot(x(i,:),Wave_U(i,:),'-b') pause(0.3) end case(2) alpha=0.05; DT=0.002;%Δt method=input("select method:Simple explicit Mehods(1)"); switch method case(1) for t=2:200 for i=1:n if(i>=2&&i<=n-1) Wave_U(t,i)=Wave_U(t-1,i)+(alpha*DT/0.05^2)*(Wave_U(t-1,i+1)-2*Wave_U(t-1,i)+Wave_U(t-1,i-1));%wave equation_Euler Explicit Method elseif(i==1) Wave_U(t,1)=Wave_U(t-1,1); else Wave_U(t,n)=Wave_U(t-1,n); end x(t,i)=x(t-1,i); end end end plot(x(1,:),Wave_U(1,:),'-r') hold on axis([-3 3 -1 2]) for i=10:20:200 plot(x(i,:),Wave_U(i,:),'-b') pause(0.3) end case(3) DT=0.02;%Δt method=input("select method:\nEuler explicit Mehods(1),Upstream Differencing Method(2),Second Order Upwind Method(3),\nLeap Frog(4),Lax-Wendroff(5),MacCormac(6)"); %选择有限差分的方法 switch method case(1)%Euler explicit Mehods for t=2:200 for i=1:n if(i>=2&&i<=n-1) Wave_U(t,i)=Wave_U(t-1,i)-(Wave_U(t-1,i)*DT/0.05)*(Wave_U(t-1,i+1)-Wave_U(t-1,i));%wave equation_Euler Explicit Method elseif(i==1) Wave_U(t,1)=Wave_U(t-1,1); else