东南大学_数值分析_第七章_偏微分方程数值解法

东南大学 数值分析 上机练习 算法与程序设计实验报告 第七章第七章 偏微分方程数值解法偏微分方程数值解法 Crank Nicolson 格式格式 学号 姓名 上机题目要求见教材 P346 10 题 一 算法原理一 算法原理 本文研究下列定解问题 抛物型方程 2 2 0 0 0 0 0 1 0 uu af x txltT tx u xxxl uttutttT MERGEFORMAT 1 的有限差分法 其中为正常数 为已知函数 且满足边界条件和初a f 始条件 关于式 1 的求解 采用离散化方法 剖分网格 构造差分格式 其中 网格剖分是将区域用两簇平行直线 0 0DxltT 0 0 i k xxihiM ttkkN 分割成矩形网格 其中分别为空间步长和时间步长 将式 1 中的 lT h MN 偏导数使用不同的差商代替 将得到不同的差分格式 如古典显格式 古典隐 格式 Crank Nicolson 格式等 其中 Crank Nicolson 格式具有更高的收敛阶数 应用更广泛 故本文采用 Crank Nicolson 格式求解抛物型方程 Crank Nicolson 格式推导 格式推导 在节点处考虑式 1 有 2 ik x t MERGEFORMAT 2 2 222 ikikik uu x tax tf x t tx 2 对偏导数用中心差分展开 2 ik u x t t 23 11 1 3 1 224 kk ikikikiikik uu x tu x tu x txtt tt MERGEFORMAT 3 第七章 偏微分方程数值解法 将在节点和表示为 2 2 2 ik u x t x ik x t 1 ik x t 222 1 222 24 1 22 1 22 8 ikikik kk iikik uuu x tx tx t xxx u xtt xt MERGEFORMAT 4 对以上两个偏导数用二阶差分展开 MERGEFORMAT 2 11 22 24 i 1 4 1 2 12 ikikikik kk ikii u x tu xtu x tu xt xh hu txx x 5 2 111111 22 24 1i 1 4 1 2 12 ikikikik kk ikii u x tu xtu x tu xt xh hu txx x MERGEFORMAT 6 将式 4 5 6 分别代入式 3 略去高阶小量 用代替并化简得 k i u ik u x t 2 1111 1111 2 1 22 22 kkkkkkkk iiiiiiiiik a uuuuuuuufx t h MERGEFORMAT 7 令 将式 7 联合式 1 初始条件和边界条件 用矩阵的形式表示为 2 rah 东南大学 数值分析 上机练习 算法与程序设计实验报告 1 11 1 22 1 22 1 11 11 22 11 2222 11 2222 11 22 kk kk kk MM kk MM rr rr uu rrrr rr uu rrrruu rr uu rr rr 11 2 2 11 22 2 2 22 kkk k Mk Mkkk r fx ttt fx t fxt r fxttt MERGEFORMAT 8 Crank Nicolson 格式的截断误差为 具有较高的精度 22 ROh 二 计算代码二 计算代码 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 三 计算结果及分析三 计算结果及分析 对于定解问题 MERGEFORMAT 9 2 2 1 0 01 01 0 01 0 1 01 x tt uu xt tx u xex uteutet 取空间步长 时间步长 采用 Crank Nicolson 格式计算 并1 xM 1 tN 将计算结果与精确值作比较 计算调用程序 Crank Nicolson m 如下 x t y xe 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 在四个目标点处的精确解 abs U CN U real err err sqrt U CN U real U CN U real 平方误差和开方 东南大学 数值分析 上机练习 算法与程序设计实验报告 M M 2 N N 2 步长反复二分 end for ii 2 kk err ii err ii 1 误差观察 end fprintf 16 15f err 在 MATLAB 中运行以上程序 计算结果如表 1 所示 表 1 Crank Nicolson 格式计算结果及其误差 40MN 输出点 x t数值解精确解误差 0 2 1 0 3 3201482661014403 3201169227365470 000031343364893 0 4 1 0 4 0552499870950874 0551999668446750 000050020250413 0 6 1 0 4 9530861218311254 9530324243951150 000053697436010 0 8 1 0 6 0496862214381036 0496474644129470 000038757025156 从表 1 可以看出 在空间步长和时间步长均为 0 25 时 采用 Crank Nicolson 格式计算抛物线方程 5 的数值解与精确解之间的误差达到了 故很好的实现偏微分方程的数值解法 且具有较高的精度和较低的运算 4 10 复杂度 同时 为了比较不同的步长所带来的计算精度 本文通过将步长反复 二分 观察四个输出点的误差均方和减小的规律 如表 2 所示 表 2 不同步长下的误差及其规律 k步长误差均方和 k err误差减少倍数 1 kk errerr 10 250 000088712758987 20 1250 0000221800894380 250021413953775 30 06250 0000055451415520 250005373866960 40 031250 0000013862901650 250000861545108 50 0156250 0000003465612770 249991874539774 60 00781250 0000000865617080 249773167584743 从表 2 可以看出 随着空间步长和时间步长的减小 数值计算的误差呈下 降趋势 说明 随着时间步长和空间步长的减小 计算结果可以收敛到精确解 同时 通过观察表 2 最后一列 可见 时间步长和空间步长同时减半 计算误 差降低为原来的四分之一 这与 Crank Nicolson 格式的截断误差 其中为时间步长 为空间步长 是符合一致的 22 ROh h