微分方程数值解第三次上机报告.doc

南京信息工程大学 实验（实习）报告 实验课程 微分方程数值解 实验名称 第三次实验 实验日期 2016 指导老师 专业 年级 姓名 学号 得分 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 实验目的： *对possion方程求解格式进行编程，并求解一个定解问题； *对迎风格式进行编程，分别适应周期边界条件的定解问题（用矩阵形式循环求解）与非周期边界条件的定解问题（元素直接循环），并各求解两个定解问题； 实验内容： 编程内容 已编程possion方程定解问题的求解程序，并用于求解实例，见possion.m文件； 编程迎风格式求解非周期边界条件的定解问题（元素直接循环），两个定解问题分别对应于，见upwind1.m与upwind2.m文件. 编程迎风格式求解周期边界条件的定解问题（用矩阵形式循环求解），两个定解问题分别对应于，见upwindm1.m与upwindm2.m文件. 问题求解 Possion方程定解问题求解： 求解如下一道possion方程定解问题： 它的解析解为： *数值求解算法： 对方向以相同步长离散化，每一行（列）均分为等份，故每一行内 点个数为个(列为个，实际上取相同步长时相等).上下左右边界条件均已知，内点求法为取上下左右四点的均值。即每个点对应一个方程，联立成方程组。再化为矩阵形式，求解方程：即如下形式： 其中 矩阵的维度都为，即行列元素个数都与每一行（列）内点个数相同。故上述A矩阵的维度为 编程时矩阵易构造，矩阵的构造可用函数实现，函数调用方法为： . 用于产生一个稀疏矩阵，其元素是矩阵中的列向量放在由行向量指定的对角线位置上（上述调用方法语句中变量名与本文问题无关）。 这里取，即下次对角线，主对角线与上次对角线三个位置。故矩阵构造程序语句如下: o=ones(N-1,1); C=-1/h^2*spdiags([o -4*o o],-1:1,N-1,N-1); 大稀疏矩阵的构造可以使用函数，其调用方法为： 意义可以理解为用矩阵的每个元素乘以矩阵，得到的矩阵替换原来元素的位置。 由矩阵构成，可分成两部分来构造矩阵，程序语句如下： p=ones(M-1,1); A1=kron(eye(M-1),C); A2=kron(spdiags([p p],[-1 1],M-1,M-1),D); A=A1+A2; 最后利用左除即可求解出所有内点。编程输出精确解三维图像，数值解三维图像和绝对误差三维分布： 图1（1）：精确解三维图 图1（2）：数值解三维图 图1（3）：绝对误差分布三维图 可以看出数值求解结果直观上很好的贴合精确解。且从绝对误差分布来看，虽然在两个分量上会有误差累积，但一直都控制在数量级。 一阶波动方程定解问题求解： 一阶波动方程定解问题一般形式如下： 随着问题中系数符号的不同会需要不同的边界条件，故求解格式上也有略微差别。只有一个边界条件时，适合对元素直接时间层逐层循环求解；当存在周期边界条件，即左右边界相等时，可以化为矩阵形式循环求解。这里求解如下四个定解问题： *非周期边界条件的的一个定解问题： （1） 其精确解为 *非周期边界条件的的一个定解问题： （2） 其精确解为 *周期边界条件的的一个定解问题： （3） 其精确解为 *周期边界条件的的一个定解问题： （4） 其精确解为 编程分别得到它们的精确解三维图，数值解三维图和绝对误差分布三维图如下： 图2（1）：问题（1）精确解 图2（2）：问题（1）数值解 图2（3）：问题（1）绝对误差分布 图3（1）：问题（2）精确解 图3（2）：问题（2）数值解 图3（3）：问题（2）绝对误差分布 图4（1）：问题（3）精确解 图4（2）：问题（3）数值解 图4（3）：问题（3）绝对误差分布 图5（1）：问题（4）精确解 图5（2）：问题（4）数值解 图5（3）：问题（4）绝对误差分布 可以看出数值解能很好地贴近精确解，直观上无法观察出差别。同时也说明算法格式稳定，这是因为四个问题的程序中，离散化方向步长均取，方