东南大学数值分析第七部分偏微分方程数值解法.docVIP

第七章 偏微分方程数值解法 ——Crank-Nicolson格式 ****(学号) *****(姓名) 上机题目要求见教材P346,10题(抛物型方程) 的有限差分法,其中为正常数,为已知函数,且满足边界条件和初始条件?关于式(1)的求解,采用离散化方法,剖分网格,构造差分格式?其中,网格剖分是将区域用两簇平行直线 分割成矩形网格,其中分别为空间步长和时间步长?将式(1)中的偏导数使用不同的差商代替,将得到不同的差分格式,如古典显格式?古典隐格式?Crank-Nicolson格式等?其中,Crank-Nicolson格式具有更高的收敛阶数,应用更广泛,故本文采用Crank-Nicolson格式求解抛物型方程? Crank-Nicolson格式推导:在节点处考虑式(1),有 对偏导数用中心差分展开 将在节点和表示为 对以上两个偏导数用二阶差分展开 将式(4)(5)(6)分别代入式(3),略去高阶小量,用代替并化简得 令,将式(7)联合式(1)初始条件和边界条件,用矩阵的形式表示为: Crank-Nicolson格式的截断误差为,具有较高的精度? 二?计算代码 Crank_Nicolson格式完整代码 function U=Crank_Nicolson(f,a,x0,xn,dx,t0,tm,dt,g,s0,sn) %Crank_Nicolson格式求解抛物线型偏微分方程 % du/dt-a*d2u/dx2=f(x,t) %Input - f 抛物方程右端函数 % - a为二阶导系数 % - x0 and xn 空间域端点 % - t0 and tm 时间域端点 % - dx为空间步长,dt为时间步长 % - g 为初始条件函数 % - s0,sn为边界条件函数 %Output - U 输出,横坐标为空间,纵坐标为时间 M=(tm-t0)/dt;N=(xn-x0)/dx;%网格数 x=x0+dx:dx:xn-dx; t=t0:dt:tm; u=feval(g,x);u=u; r=a*dt/dx/dx;%步长比 %Crank-Nicolson格式:A*u_(k+1)=B*u_k+C A=-r/2*[zeros(1,N-1);eye(N-2,N-2),zeros(N-2,1)]-r/2* ... [zeros(N-2,1),eye(N-2,N-2);zeros(1,N-1)]+(1+r)*eye(N-1,N-1); A=inv(A); B=r/2*[zeros(1,N-1);eye(N-2,N-2),zeros(N-2,1)]+r/2* ... [zeros(N-2,1),eye(N-2,N-2);zeros(1,N-1)]+(1-r)*eye(N-1,N-1); U=[]; for k=1:M C=dt*feval(f,x,t(k)+0.5*dt); C=C; C(1)=C(1)+r/2*(feval(s0,t(k))+feval(s0,t(k+1))); C(N-1)=C(N-1)+r/2*(feval(sn,t(k))+feval(sn,t(k+1))); u=A*(B*u+C); U=[U;u]; end 三?计算结果及分析 对于定解问题 取空间步长,时间步长,采用Crank_Nicolson格式计算,并将计算结果与精确值作比较,计算调用程序Crank_Nicolson.m,如下 clc format long f=inline(0*x+0*t);%抛物方程右端函数 g=inline(exp(x));%初始条件函数 s0=inline(exp(t));sn=inline(exp(t+1));%边界条件函数 a=1;x0=0;xn=1;t0=0;tm=1; M=40;N=40; err=[]; kk=4; for ii=1:kk dx=1/N;dt=(tm-t0)/M;%空间步长,时间步长 U=Crank_Nicolson(f,a,x0,xn,dx,t0,tm,dt,g,s0,sn);%调用Crank_Nicolson函数 U_CN=U(length(U),N/5:N/5:4*N/5); %在四个目标点处的数值解 fu=inline(exp(x+t)); %精确方程 xx=0.2:0.2:0.8;tt=1.0; U_real=feval(fu,xx,tt); %在四个