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

第七章 偏微分方程数值解法 东南大学《数值分析》上机练习——算法与程序设计实验报告 第七章 偏微分方程数值解法 ——Crank-Nicolson格式 ****（学号） *****（姓名） 上机题目要求见教材P346,10题。 一、算法原理 本文研究下列定解问题(抛物型方程)  MACROBUTTON MTPlaceRef \* MERGEFORMAT  SEQ MTEqn \h \* MERGEFORMAT ( SEQ MTEqn \c \* Arabic \* MERGEFORMAT 1) 的有限差分法，其中为正常数，为已知函数，且满足边界条件和初始条件。关于式(1)的求解，采用离散化方法，剖分网格，构造差分格式。其中，网格剖分是将区域用两簇平行直线 分割成矩形网格，其中分别为空间步长和时间步长。将式(1)中的偏导数使用不同的差商代替，将得到不同的差分格式，如古典显格式、古典隐格式、Crank-Nicolson格式等。其中，Crank-Nicolson格式具有更高的收敛阶数，应用更广泛，故本文采用Crank-Nicolson格??求解抛物型方程。 Crank-Nicolson格式推导：在节点处考虑式(1)，有  MACROBUTTON MTPlaceRef \* MERGEFORMAT  SEQ MTEqn \h \* MERGEFORMAT ( SEQ MTEqn \c \* Arabic \* MERGEFORMAT 2) 对偏导数用中心差分展开  MACROBUTTON MTPlaceRef \* MERGEFORMAT  SEQ MTEqn \h \* MERGEFORMAT ( SEQ MTEqn \c \* Arabic \* MERGEFORMAT 3) 将在节点和表示为  MACROBUTTON MTPlaceRef \* MERGEFORMAT  SEQ MTEqn \h \* MERGEFORMAT ( SEQ MTEqn \c \* Arabic \* MERGEFORMAT 4) 对以上两个偏导数用二阶差分展开  MACROBUTTON MTPlaceRef \* MERGEFORMAT  SEQ MTEqn \h \* MERGEFORMAT ( SEQ MTEqn \c \* Arabic \* MERGEFORMAT 5)  MACROBUTTON MTPlaceRef \* MERGEFORMAT  SEQ MTEqn \h \* MERGEFORMAT ( SEQ MTEqn \c \* Arabic \* MERGEFORMAT 6) 将式(4)(5)(6)分别代入式(3)，略去高阶小量，用代替并化简得  MACROBUTTON MTPlaceRef \* MERGEFORMAT  SEQ MTEqn \h \* MERGEFORMAT ( SEQ MTEqn \c \* Arabic \* MERGEFORMAT 7) 令，将式(7)联合式(1)初始条件和边界条件，用矩阵的形式表示为：  MACROBUTTON MTPlaceRef \* MERGEFORMAT  SEQ MTEqn \h \* MERGEFORMAT ( SEQ MTEqn \c \* Arabic \* MERGEFORMAT 8) 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-Nicolso