一维扩散方程的有限差分法matlab

1一维扩散方程的有限差分法——计算物理实验作业七陈万物理学2013级13020011006题目：编程求解一维扩散方程的解)()0(|),()0,0(02221110max022axcnubuaxcnubuaetxuttaxxuDtuxt取1.0,1.0,1.0,10,0.1,0,1,1,0,1,1max0222111hDtacbacba。输出t=1,2,...,10时刻的x和u(x),并与解析解u=exp(x+0.1t)作比较。主程序：%一维扩散方程的有限差分法clear,clc;%定义初始常量a1=1;b1=1;c1=0;a2=1;b2=-1;c2=0;a0=1.0;t_max=10;D=0.1;h=0.1;tao=0.1;%调用扩散方程子函数求解u=diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2);子程序1：functionoutput=diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2)%一维扩散方程的有限差分法，采用隐式六点差分格式(Crank-Nicolson)%a0:x的最大值%t:_max:t的最大值%h:空间步长%tao:时间步长%D：扩散系数%a1,b1,c1是（x=0）边界条件的系数；a2,b2,c2是（x=a0）边界条件的系数x=0:h:a0;2n=length(x);t=0:tao:t_max;k=length(t);P=tao*D/h^2;P1=1/P+1;P2=1/P-1;u=zeros(k,n);%初始条件u(1,:)=exp(x);%求A矩阵的对角元素dd=zeros(1,n);d(1,1)=b1*P1+h*a1;d(2:(n-1),1)=2*P1;d(n,1)=b2*P1+h*a2;%求A矩阵的对角元素下面一行元素ee=-ones(1,n-1);e(1,n-1)=-b2;%求A矩阵的对角元素上面一行元素ff=-ones(1,n-1);f(1,1)=-b1;R=zeros(k,n);%求R%追赶法求解fori=2:kR(i,1)=(b1*P2-h*a1)*u(i-1,1)+b1*u(i-1,2)+2*h*c1;forj=2:n-1R(i,j)=u(i-1,j-1)+2*P2*u(i-1,j)+u(i-1,j+1);endR(i,n)=b2*u(i-1,n-1)+(b2*P2-h*a2)*u(i-1,n)+2*h*c2;M=chase(e,d,f,R(i,:));u(i,:)=M';plot(x,u(i,:));axis([0a00t_max]);pause(0.1)endoutput=u;%绘图比较解析解和有限差分解[X,T]=meshgrid(x,t);Z=exp(X+0.1*T);surf(X,T,Z),xlabel('x'),ylabel('t'),zlabel('u'),title('解析解');figure3surf(X,T,u),xlabel('x'),ylabel('t'),zlabel('u'),title('有限差分解');子程序2：functionM=chase(a,b,c,f)%追赶法求解三对角矩阵方程，Ax=f%a是对角线下边一行的元素%b是对角线元素%c是对角线上边一行的元素%M是求得的结果，以列向量形式保存n=length(b);beta=ones(1,n-1);y=ones(1,n);M=ones(n,1);fori=(n-1):(-1):1a(i+1)=a(i);end%将a矩阵和n对应beta(1)=c(1)/b(1);fori=2:(n-1)beta(i)=c(i)/(b(i)-a(i)*beta(i-1));endy(1)=f(1)/b(1);fori=2:ny(i)=(f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1));endM(n)=y(n);fori=(n-1):(-1):1M(i)=y(i)-beta(i)*M(i+1);endend结果：4对比分析两图，结果令人满意。取出t_max时刻的u值分析：（111~uu）有限差分解如下：2.7182818283.0041660243.3201169233.6692966684.0551999674.4816890704.9530324245.4739473926.0496474646.6858944427.389056099解析解如下：2.7145201222.9995109753.3144223293.6623954644.0469014544.4717757924.9412566525.460027166.0332621166.6666796077.36659805分析数据可知误差量级为)(O)(O22h=0.12+0.12=0.02总结：（1）隐式六点差分格式(Crank-Nicolson)基本思想是用前一时刻的三个点表示后一时刻的三个点。因为不是直接表示u(k+1,i)，故称为隐式差分。（2）x由等步长被分割为N个点，列出N个方程。采用追赶法求解得到结果。原理很简单，关键是求解AU=R中的A和R。（3）子函数2的功能是实现追赶法，该程序中没有直接用A来表示三对角矩阵，而是把3列对角元素直接拿出来，因此在调用时应当注意各对角元素的位置，避免调用错误。