偏微分实验报告1

偏微分方程数值解法第三次上机实验报告一、实验题目：用线性元求解下列问题的数值解：2,1,1(x,1)u(x,1)0,1x 1(1,y)1,u (1,y)0,11xx u x y u u y =--<<⎧⎪-==-<<⎨⎪-==-<<⎩ （精确到小数点后四位）二、实验过程：利用PDEToolbox 工具箱求解该偏微分方程。

分析：方程是Possion 方程形式c u au f -+=，其中c=-1，a=0，f=-2； 第一个边界条件是Dirichlet 条件，第二个边界条件是Neumann 条件。

1.在MA TLAB 命令窗口键入pdetool 并运行，打开PDEToolbox 界面；2.在Options 菜单下选择Grid 命令，显示网格，能更容易确定所绘图形的大小；3.绘出区域，选择Boundary 的Boundary Mode ，双击每个边界，设置边界相应的参数值；4.选择PDE 菜单中PDE Mode 命令，进入PDE 模式。

单击PDE 菜单中PDE Specification ….选项，设置方程类型及参数；5.选择Mesh 菜单中Initialize Mesh 命令，进行网格剖分，再选择Refine Mesh 命令，进行网格加密，如下图：三、实验结果：选择Solve 菜单中solve PDE 命令，解偏微分方程，其图形解如图：图1 图形解图2 三维图形解图3 解的等值线图和矢量场图选择Mesh菜单中的Export Mesh，得到结点xy坐标；选择Solve菜单中的Export Solution…，得到每个节点处的值，输出u，即解的数值。

四、实验总结：通过本次试验，掌握了利用Matlab中的PDE求解工具得到PDE的解的方法，并对偏微分方程的形式有了更多的掌握。

偏微分方程报告范文偏微分方程是研究多变量函数的微分方程，其中函数的未知量既依赖于自变量，又依赖于多个自变量。

偏微分方程在物理学、工程学和经济学等领域中有着广泛的应用。

本报告将介绍偏微分方程的基本概念、求解方法和应用领域。

一、偏微分方程的基本概念偏微分方程是由未知函数的偏导数和自变量构成的方程。

常见的偏微分方程包括波动方程、传热方程和扩散方程等。

偏微分方程根据阶数可分为一阶和二阶偏微分方程。

一阶偏微分方程中只涉及到未知函数的一阶偏导数，一般可以通过变量分离的方法求解。

二阶偏微分方程中涉及到未知函数的二阶偏导数，求解方法一般包括分离变量法、特征线法和变换法等。

二、偏微分方程的求解方法1.分离变量法：假设未知函数可以表示为两个只依赖于单个自变量的函数的乘积形式，然后将该形式代入到偏微分方程中，再将方程两边关于不同的自变量求积分，从而得到方程的通解。

2.特征线法：通过特征线曲线的方法将偏微分方程转化为常微分方程。

先找出特征线曲线，然后在特征线上引入新的变量，使得偏微分方程变为常微分方程，进而求解。

3.变换法：通过适当的变量变换，将原偏微分方程转化为一个更容易求解的形式。

常用的变换方法有坐标变换、函数变换和变量替换等。

三、偏微分方程的应用领域1.物理学：偏微分方程在物理学中有着广泛的应用。

例如，波动方程可以描述声波、光波和电磁波等在介质中的传播；传热方程可以描述热传导过程；薛定谔方程和波恩-奥本海默方程可以描述量子力学中的粒子行为等。

2.工程学：偏微分方程在工程学中被广泛应用于流体力学、结构力学和电磁场等领域。

例如，纳维-斯托克斯方程用于描述流体的运动；弹性方程用于描述结构的变形和应力分布等。

3.经济学：偏微分方程在经济学中应用较多，尤其是在金融学中。

例如，布莱克-斯科尔斯方程用于定价期权；黑-舒尔斯方程用于描述衍生品的定价和风险管理等。

通过对偏微分方程的研究和求解，可以更好地理解自然界的现象和规律，并为解决实际问题提供数学模型和解决方法。

预备实验偏微分方程工具箱的功能和利用一、实验目的PDE Toolbox提供了研究和求解空间二维偏微分方程问题的一个壮大而又灵活有效的环境。

本次实验了解PDE Toolbox的功能和PDE图形用户界面。

PDE Toolbox的功能包括：设置PDE 定解问题，即设置二维定解区域、边界条件和方程的形式和系数；用有限元法求解PDE，即网格的生成、方程的离散和求出数值解；解得可视化。

二、实验内容一、解偏微分方程的一个例子解Poisson方程：f∆-，边界条件为齐次dirichlet类型。

u=第一步：启动MATLAB，键入pdetool，按回车键确信即可启动GUI(图形用户界面)，然后在options菜单下选择grid命令，打开栅格。

栅格的利用，能利用户容易确信所画图形的大小。

第二步：散布完成平面几何造型：R1-C1-E1+R2+C2。

用菜单或快捷工具别离画矩形R1，矩形R2，椭圆E1，圆C1，圆C2。

画圆时，第一选中椭圆工具，按鼠标右键并拖动即可，或在按Ctrl的同时，拖动鼠标也可绘制图。

然后在set formula栏，进行编辑并用算术运算符将图形对象名称连接起来，或删除默许的表达式直接键入R1-C1-E1+R2+C2。

选择Boundary菜单中Boundary Mode 命令，进入边界模式。

单击Boundary 菜单中Remove All Subdomain Borders选项，去掉子域边界，若是想将几何信息和边界信息进行存储，应选择Boundary菜单中的Export Decomposed Geometry, Boundary Cond’s…命令，将它们别离贮存于g,b变量中，通过MATLAB形成M文件。

第三步：选取边界，单击Boundary菜单中Specify Boundary Conditions…选项，打开Boundary Conditions对话框，输入边界条件，本例取缺省条件，即将全数边界设为齐次dirichlet条件，边界颜色显示为红色。

微分方程数值解实验指导（一）实验一: 二阶椭圆型方程差分格式的程序实现1. 实验内容用五点差分格式求解Poisson 方程边值问题⎩⎨⎧∂∈=∈=-=∆,),(,0,),(,:2G y x u G y x f u （1） 其中}1|||,||),{(<=y x y x D 。

（1）用正方形网格)(h k =列出相应的差分方程。

（2）对161,81,51,21=h 分别求出数值解，观察数值解的情况，分析计算结果。

（最好画出数值解的图形）注意：在区域G 的边界上为齐次Dirichlet 条件，在这类边界上不需要给出差分格式。

2. 实验目的及要求按照给定的差分格式编程实现求出数值解；结合格式的相容性和收敛性条件简单分析计算结果。

要求在实验课上算出数值结果；按要求格式写出实验报告；下次实验课前交本次实验的实验报告。

3. 实验原理与实验过程: 以下是求解问题的步骤第一步: 对求解区域作网格剖分。

按照要求得网格剖分图如下（图1为21=h 的情况，图2为51=h 的情况）-1-0.8-0.6-0.4-0.20.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.81区域划分示意图图1-1-0.8-0.6-0.4-0.200.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.81区域划分示意图图2第二步: 差分法的目的是：求边值问题的解),(y x u 在节点),(j i y x 的近似值ij u （m j n i ,...,2,1,,2,1==， ）。

为此，需要构造逼近微分方程定解问题的差分格式。

采用五点差分格式，取定x 轴和y 轴方向的步长相等，即h=k ，作两族与坐标轴平行的直线：ih x i +-=1，i=0,1,…,2/h,jh y j +-=1，j=0,1,…,2/h,两族直线的交点为网点（或节点）。

位于G 内的节点为内点，位于G 的边界上的节点为界点。

偏微分方程数值实验报告八实验题目：利用有限差分法求解.0)1(,0)1(),()()(==-=+''-u u x f x u x u 真解为)1()(22x e x u x -=-实现算法：对于两点边值问题,)(,)(,,d 22βα==∈=-b u a u l x f dxu(1)其中),(b a l =f b a ),(<为],[b a l =上的连续函数，βα，为给定常数.其相应的有限差分法的算法如下：1.对求解区域做网格剖分，得到计算网格.在这里我们对区间l 均匀剖分n 段，每个剖分单元的剖分步长记为nab h -=.2.对微分方程中的各阶导数进行差分离散，得到差分方程.运用的离散方法有：方法一:用待定系数和泰勒展开进行离散)()()()(d )(d 111122++--++≈i i i i i i i i x u x u x u x x u ααα方法二：利用差商逼近导数21122)()(2)()(d )(d h x u x u x u x x u i i i i i -++-≈(2)将(2)带入(1)可以得到)()()()(2)(211u R x f h x u x u x u i i i i i +=+---+,其中)(u R i 为无穷小量，这时我们丢弃)(u R i ，则有在i x 处满足的计算公式：1,...,1)()()(2)(211-==+---+n i x f hx u x u x u i i i i ，(3)3.根据边界条件，进行边界处理.由(1)可得βα==n u u ,0(4)称(3)(4)为逼近(1)的差分方程，并称相应的数值解向量1-n U 为差分解，i u 为)(i x u 的近似值.4.最后求解线性代数方程组，得到数值解向量1-n U .实验题目：用Lax-Wendroff 格式求解方程：.4sin 1),1(],1,0[,2sin 1)0,(,0),1,0(,02t t u x x x u t x u u x t ππ+=∈+=>∈=- （1） （精确解).2(2sin 1t x u ++=π） 数值边值条件分别为： （a ）).(20101n 0nn nu u hu u -+=+τ （b ）.1n 0nu u =（c ）.02-12111n 0=++++n n u u u请将计算结果与精确解进行比较。

偏微分方程数值解法上机报告（一）一、实验题目：用Ritz-Galerkin 方法求解边值问题2u '',01(0)0,(1)1u x x u u ⎧-+=<<⎨==⎩的第n 次近似()n u x ，基函数()sin(),1,2,...,i x i x i n ϕπ==.二、实验目的：通过本次上机实验，理解求解初值问题的变分问题的最重要的近似解法——Ritz-Galerkin 方法，以便为学习有限元法打好基础。

此外，要熟悉用Matlab 解决数学问题的基本编程方法，提高运用计算机解决问题的能力。

三、实验代码：n=5;syms x;for i=1:np(i)=sin(i*pi*x);q(i)=-i^2*pi^2*sin(i*pi*x);endfor i=1:nb(i)=2*int(p(i),0,1);for j=1:nA(i,j)=int((-q(j)+p(j))*p(i),0,1);endendt=inv(A)*b'四、运行结果：t=2251799813685248/3059521645650671/pi281474976710656/9481460623939047/pi281474976710656/43582901062631895/pi五、总结：通过本次上机，我了解了Ritz-Galerkin 方程 n j j p f cj p i p a n i i ,...,2,1)),(,())(),((1==∑=，明白了用Ritz-Galerkin 方法解决边值问题的变分问题的基本原理，并接近一步提高自己的编程动手能力，受益匪浅。

偏微分方程数值解法上机报告（二）一、 实验题目：用线性元求下列边值问题的数值解2''2sin ,0142(0)0,'(1)0y y x x y y ππ⎧-+=<<⎪⎨⎪==⎩二、 实验目的：通过本次上机，熟悉和掌握用Galerkin 法观点出发导出的求解处置问题数值解的线性有限元法。

微分方程数值解实验指导（一）实验一: 二阶椭圆型方程差分格式的程序实现1. 实验内容用五点差分格式求解Poisson 方程边值问题⎩⎨⎧∂∈=∈=-=∆,),(,0,),(,:2G y x u G y x f u （1） 其中}1|||,||),{(<=y x y x D 。

（1）用正方形网格)(h k =列出相应的差分方程。

（2）对161,81,51,21=h 分别求出数值解，观察数值解的情况，分析计算结果。

（最好画出数值解的图形）注意：在区域G 的边界上为齐次Dirichlet 条件，在这类边界上不需要给出差分格式。

2. 实验目的及要求按照给定的差分格式编程实现求出数值解；结合格式的相容性和收敛性条件简单分析计算结果。

要求在实验课上算出数值结果；按要求格式写出实验报告；下次实验课前交本次实验的实验报告。

3. 实验原理与实验过程: 以下是求解问题的步骤第一步: 对求解区域作网格剖分。

按照要求得网格剖分图如下（图1为21=h 的情况，图2为51=h 的情况）-1-0.8-0.6-0.4-0.20.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.81区域划分示意图图1-1-0.8-0.6-0.4-0.200.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.81区域划分示意图图2第二步: 差分法的目的是：求边值问题的解),(y x u 在节点),(j i y x 的近似值ij u （m j n i ,...,2,1,,2,1==， ）。

为此，需要构造逼近微分方程定解问题的差分格式。

采用五点差分格式，取定x 轴和y 轴方向的步长相等，即h=k ，作两族与坐标轴平行的直线：ih x i +-=1，i=0,1,…,2/h,jh y j +-=1，j=0,1,…,2/h,两族直线的交点为网点（或节点）。

位于G 内的节点为内点，位于G 的边界上的节点为界点。

偏微分中心差分格式实验报告实验目的：1.掌握偏微分的中心差分格式；2.理解中心差分格式的精度和稳定性。

实验原理：中心差分是一种常用的数值求解偏微分方程的格式，其基本思想是用函数在两个点的导数的平均值来近似函数在这两个点中间的导数值。

具体来说，对于一维的偏微分方程，中心差分格式可以表述为：f'(x)=(f(x+h)-f(x-h))/(2h)其中f'(x)表示x点处的导数，h表示步长。

实验步骤：1.编写一个计算函数在任意给定点x处的导数值的中心差分程序；2.给定一个函数f(x)，例如f(x)=x^2，计算在一定范围内的该函数在每个点处的导数值；3.比较计算的导数值与理论值的差异，并分析中心差分格式的精度；4.对给定步长h，逐渐减小h，计算导数值，并观察数值的变化，分析中心差分格式的稳定性。

实验结果与分析：以函数f(x)=x^2为例，给定步长h=0.1，计算在范围[-1,1]内的函数f(x)在每个点处的导数值。

实验结果如下表所示：x，f'(x)，理论值，误差-1.0，-1.999，-2，0.001 -0.9，-1.899，-1.8，0.099 -0.8，-1.698，-1.6，0.098 -0.7，-1.397，-1.4，0.003 -0.6，-0.996，-1，0.004 -0.5，-0.495，-0.5，0.005 -0.4，0.204，0，0.204-0.3，0.615，0.6，0.015 -0.2，1.216，1.2，0.016 -0.1，1.797，1.8，0.003 0.0，1.996，2，0.0040.1，2.193，2.2，0.007 0.2，2.792，2.8，0.008 0.3，3.293，3.3，0.007 0.4，3.594，3.6，0.006 0.5，3.896，3.9，0.004 0.6，4.437，4.4，0.037 0.7，4.998，5，0.0020.9，6.795，6.8，0.0051.0，7.993，8，0.007从实验结果可以看出，随着x的增大，计算的导数值与理论值之间的误差也在增大，但整体上相对较小。

二阶常微分方程的中心差分求解学校:中国石油大学（华东）理学院 姓名：张道德一、 实验目的1、 构造二阶常微分边值问题:22,(),(),d uLu qu f a x bdx u a u b αβ⎧=-+=<<⎪⎨⎪==⎩其中,q f 为[,]a b 上的连续函数，0;,q αβ≥为给定常数的中心差分格式；2、 根据中心差分格式求解出特定例题的数值解，并与该例题的精确解进行比较。

二、 中心差分格式的构造将区间[,]a b 分成N 等分，分点为： 0,1,2,,i x a ih i N =+=()/h b a N =-。

于是我们得到区间的一个网络剖分。

称为网格的节点称为步长。

得中心差分格式为：11202,1,2,,1,,.i i i h i i i i N u u u L u q u f i N h u u αβ+--+⎧=-+==-⎪⎨⎪==⎩其中式中(),()i i i i q q x f f x ==。

三、 差分格式求解根据中心差分格式可以构造出：1112222222233322212211210012101201001200N N N u f q h h u f q h h h u f q h hh q u f h h ---⎡⎤⎡⎤⎡⎤+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=-+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦可以看出系数矩阵为三对角矩阵，而对于系数矩阵为三对角矩阵的方程组可以用“追赶法”求解，则可以得出二阶常微分方程问题的数值解。

四、 举例求解我们选取的二阶常微分方程边值问题为：222242,01(0)1,(1),x d u Lu x u e x dx u u e ⎧=-+=-<<⎪⎨⎪==⎩其精确解为：2x u e=。

则我们具体求解出的解如下：1、 当10N =时，数值解与精确解为： (1) 表1、10N =时，数值解与精确解统计表x 的值 0.10.20.30.40.5 u 的数值解 1.011069 1.042744 1.096904 1.176896 1.28789 u 的精确解 1.01005 1.040811 1.094174 1.173511 1.284025 两者之差 0.001019 0.001934 0.002729 0.003385 0.003864 x 的值 0.60.70.80.9u 的数值解 1.437443 1.636363 1.900001 2.250209 u 的精确解 1.433329 1.632316 1.896481 2.247908 两者之差0.0041140.0040460.003520.002301将两者绘于同一图中如下：(2)结论：可以看出数值解与精确解之间的误差很小， 在 210-这样一个数量级上。

图像处理偏微分方法实验报告偏微分方程在图像处理领域的研究开始于六、七十年代，从最初的去噪角度和图像恢复的角度出发，相应地引入了偏微分方程。

但直到九十年代才比较系统地将偏微分方程引入图像处理领域，结合其他一些数学工具如数学形态学和仿射几何等，形成了比较完整的理论体系。

由s．Osher等提出的水平集方法具有良好的几何插值性，将灰度插值转化为曲线插值，在图像处理中有很大的影响，在图像修复与去噪、边缘检测、图像匹配、图像识别等方面都取得了相对较好的结果。

一、偏微分方程在图像处理中的应用偏微分方程在图像处理领域的研究开始于六、七十年代，Gabor，Rudin，Osher 等人分别从去噪的角度出发，把偏微分方程引入了图像处理领域。

同时，Osher 还从图像恢复的角度出发，使用最小化全局变差(Total Vadation)方法和变分法，也使用了相应的偏微分方程。

但直到上世纪九十年代，Alvarez，Lions，Morel等才首次系统地将偏微分方程引入了图像处理领域，在理论上有比较大的突破，并结合其他数学工具，如数学形态学、变分法、逼近论、仿射几何等，形成了比较完整的理论体系。

由于理论体系本身的优点，从二维图像(静止图像，如照片)到三维图像(运动图像，如电影)应用的拓展相对容易一些。

这也是偏微分方程在图像处理研究领域中的一个优势。

当然，时间这一维数有它本身的特点，与图像平面有本质的区别，在应用中需要加以注意。

概括地说，偏微分方程的图像处理技术具有以下的特点：(1)基于PDE的方法具有良好的数学基础，可以提供深刻的理论结果，并且算法具备良好的稳定性；(2)一些经典的方法如高斯滤波、中值滤波、膨胀和腐蚀在PDE的统一框架下得到全新的解释；(3)这种视角产生了新的方法，它们比经典的方法包容更多的不变性，如保持结构的滤波、线形增强等；(4)PDE是连续的模型，与具体的离散网格无关，并且具有旋转不变性。

偏微分方程在图像处理中的研究主要集中在图像处理的两个基础部分：图像滤波和图像边缘提取。

2009级数学与应用数学和信息与计算科学专业偏微分方程数值解上机实验实验题目利用有限元方法和有限差分方法求解偏微分方程完成日期 2012年12月17日学生姓名张灵刚所在班级 1102090 任课教师王晓东西北工业大学理学院应用数学系目录一．实验目的 (2)二．实验要求 (2)三．实验题目 (3)四．实验二 (4)1.实验内容 (4)2.实验原理.................................................（4）. 3算法流程. (5)4结果分析 (5)5总结讨论 (6)6源程序 (6)五. 实验三1.实验内容 (17)2.实验原理...............................................（17）. 3算法流程. (18)4结果分析……………………………………….….（18）.5总结讨论 (21)6源程序 (21)偏微分方程数值解上机实验报告实验地点：数学系机房实验时间：第13—15周，周一、四下午5、6节实验分数：占期末考试成绩的30%一、实验目的及意义掌握有限元方法和有限差分方法的程序实现；学会选择合适的有限差分格式求解一维非线性对流占优的非定常对流扩散问题；学会使用三角线性元和四边形线性元的有限元方法求解二维椭圆方程边值问题，并对计算结果进行收敛性分析；尝试采用有限元方法或有限差分方法实现二维初边值抛物型方程的大规模数值求解。

通过实验可以提高学生的动手能力，加深学生对算法的理解。

二、实验要求在下列给出的三个问题中，最少选择两个问题进行编程实现。

要求给出格式的推导过程、算法流程、实现程序、选取的网格参数、以表格或图形的方式给出计算结果、对计算结果进行分析、最后对实验进行总结和讨论。

问题2：用三角线性元和四边形线性元的有限元方法求解方程：28cos(2)sin(2),(,)(0,1)(0,1)(,0)(,1)0;(0,)(1,)sin(2).u x y x y u x u x u y u y y ππππ-∆=∈⨯====取=1/4, 1/8, 1/16, 1/32, 1/64,h 比较两种方法的计算精度，并给出数值收敛率.问题3：选用合适的数值方法求解方程：22122(444),(,)(0,1)(0,1)(0,,)(1,,)(,0,)(,1,)0;(,,0)sin()cos().y y ux y u x y tu y t u y t u x t u x t u x y x y ππππ-∂=-++∆∈⨯∂''=====求0.1,0.5,1.0t =时，点3331(,)6464、1517(,)6464、4749(,)6464、7137(,)128128、4367(,)128128、10989(,)128128、127129(,)256256、6391(,)256256、3133(,)256256处的数值解。

偏微分方程数值解上机实验报告（一）实验一一、上机题目：用线性元求解下列边值问题的数值解：-y′′+y24y=y22sin y2y，0<x<1y(0)=0，y′(1)=0二、实验程序：function S=bzx=fzero(@zfun,1);[t y]=ode45(@odefun,[0 1],[0 x]); S.t=t;S.y=y;plot(t,y)xlabel('x:´从0一直到1')ylabel('y')title('线性元求解边值问题的数值解') function dy=odefun(x,y)dy=[0 0]';dy(1)=y(2);dy(2)=(pi^2)/4*y(1)-((pi^2)/2)*sin(x*pi/2);function z=zfun(x);[t y]=ode45(@odefun,[0 1],[0 x]);z=y(end)-0;三、实验结果：1．以步长h=0.05进行逐步运算，运行上面matlab程序结果如下：2．在0<x<1区间上，随着x的不断变化，x，y之间关系如下图所示：（二）实验二四、 上机题目：求解Helmholtz 方程的边值问题：21u k u -∆-=，于(0,1)*(0,1)Ω=0u =，于1{0,01}{01,1}x y x y Γ==≤≤≤≤=U 12{0,01}{01,1}0,{01,0}{1,01}x y x y u x y x y n Γ==≤≤≤≤=∂=Γ=≤≤==≤≤∂U U 于 其中k=1,5,10,15,20五、实验程序：（采用有限元方法，这里对[0,1]*[0,1]采用n*n的划分，n为偶数）n=10;a=zeros(n);f=zeros(n);b=zeros(1,n);U=zeros(n,1);u=zeros(n,1);for i=2:na(i-1,i-1)=pi^2/(12*n)+n;a(i-1,i)= pi^2/(24*n)-n;a(i,i-1)= pi^2/(24*n)-n;for j=1:nif j==i-1a(i,j)=a(i,i-1);else if j==ia(i-1,j-1)=2*a(i-1,i-1);else if j==i+1a(i,j)=a(i,i+1);elsea(i,j)=0;endendendendenda(n,n)=pi^2/(12*n)+n;for i=2:nf(i-1,i)=4/pi*cos((i-1)*pi/2/n)-8*n/(pi^2)*sin(i*pi/2/n)+8*n/(pi^2)*sin((i-1)*pi/2/n); endfor i=1:nf(i,i)=-4/pi*cos(i*pi/2/n)+8*n/(pi^2)*sin(i*pi/2/n)-8*n/(pi^2)*sin((i-1)*pi/2/n); end%b(j)=f(i-1,j)+f(i,j)for i=1:(n-1)b(i)=f(i,i)+f(i,i+1);endb(n)=f(n,n);tic;n=20;can=20;s=zeros(n^2,10);h=1/n;st=1/(2*n^2);A=zeros((n+1)^2,(n+1)^2);syms x y;for k=1:1:2*n^2s(k,1)=k;q=fix(k/(2*n));r=mod(k,(2*n));if (r~=0)r=r;else r=2*n;q=q-1;endif (r<=n)s(k,2)=q*(n+1)+r;s(k,3)=q*(n+1)+r+1;s(k,4)=(q+1)*(n+1)+r+1;s(k,5)=(r-1)*h;s(k,6)=q*h;s(k,7)=r*h;s(k,8)=q*h;s(k,9)=r*h;s(k,10)=(q+1)*h;elses(k,2)=q*(n+1)+r-n;s(k,3)=(q+1)*(n+1)+r-n+1;s(k,4)=(q+1)*(n+1)+r-n;s(k,5)=(r-n-1)*h;s(k,6)=q*h;s(k,7)=(r-n)*h;s(k,8)=(q+1)*h;s(k,9)=(r-n-1)*h;s(k,10)=(q+1)*h;endendd=zeros(3,3);B=zeros((n+1)^2,1);b=zeros(3,1);for k=1:1:2*n^2L(1)=(1/(2*st))*((s(k,7)*s(k,10)-s(k,9)*s(k,8))+(s(k,8)-s(k,10))*x+(s(k,9)-s(k,7))*y);L(2)=(1/(2*st))*((s(k,9)*s(k,6)-s(k,5)*s(k,10))+(s(k,10)-s(k,6))*x+(s(k,5)-s(k,9))*y);L(3)=(1/(2*st))*((s(k,5)*s(k,8)-s(k,7)*s(k,6))+(s(k,6)-s(k,8))*x+(s(k,7)-s(k,5))*y);for i=1:1:3for j=i:3d(i,j)=int(int(((((diff(L(i),x))*(diff(L(j),x)))+((diff(L(i),y))*(diff(L(j),y))))-((can^2)*L(i)*L(j))),x,0, 1),y,0,1);d(j,i)=d(i,j);endendfor i=1:1:3for j=1:1:3A(s(k,(i+1)),s(k,(j+1)))=A(s(k,(i+1)),s(k,(j+1)))+d(i,j);endendfor i=1:1:3b(i)=int(int((L(i)),x,0,1),y,0,1);B(s(k,(i+1)),1)=B(s(k,(i+1)),1)+b(i);endendM=zeros((n+1)^2,n^2);j=n^2;for i=(n^2+n):-1:1if ((mod(i,(n+1)))~=1)M(:,j)=A(:,i);j=j-1;else continueendendpreanswer=M\B;answer=zeros((n+1)^2,1);j=1;for i=1:1:(n^2+n)if ((mod(i,(n+1)))~=1)answer(i)=preanswer(j);j=j+1;else answer(i)=0;endendZ=zeros((n+1),(n+1));for i=1:1:(n+1)^2s=fix(i/(n+1))+1;r=mod(i,(n+1));if(r==0)r=n+1;s=s-1;elseendZ(r,s)=answer(i);end[X,Y]=meshgrid(1:-h:0,0:h:1);surf(X,Y,Z);toc;t=toc;K=a ;B=b';U=inv(K)*Bfor i=1:nu(i,1)=4/(pi^2)*sin(pi*i/n/2);endue=U-u六、实验结果：程序中的变量can为问题中的k，为了方便比较，采用了画图的方式。

偏微分方程数值解上机实验报告（一）实验一一、上机题目：用线性元求解下列边值问题的数值解：-y′′+π24y=π22sinπ2x，0<x<1y(0)=0，y′(1)=0二、实验程序：function S=bzx=fzero(@zfun,1);[t y]=ode45(@odefun,[0 1],[0 x]);S.t=t;S.y=y;plot(t,y)xlabel('x:´从0一直到1')ylabel('y')title('线性元求解边值问题的数值解')function dy=odefun(x,y)dy=[0 0]';dy(1)=y(2);dy(2)=(pi^2)/4*y(1)-((pi^2)/2)*sin(x*pi/2);function z=zfun(x);[t y]=ode45(@odefun,[0 1],[0 x]);z=y(end)-0;三、实验结果：1．以步长h=0.05进行逐步运算，运行上面matlab程序结果如下：2．在0<x<1区间上，随着x 的不断变化，x ，y 之间关系如下图所示：（二）实验二四、 上机题目：求解Helmholtz 方程的边值问题：21u k u -∆-=，于(0,1)*(0,1)Ω=0u =，于1{0,01}{01,1}x y x y Γ==≤≤≤≤= 12{0,01}{01,1}0,{01,0}{1,01}x y x y u x y x y n Γ==≤≤≤≤=∂=Γ=≤≤==≤≤∂于其中k=1,5,10,15,20五、实验程序：（采用有限元方法，这里对[0,1]*[0,1]采用n*n的划分，n为偶数）n=10;a=zeros(n);f=zeros(n);b=zeros(1,n);U=zeros(n,1);u=zeros(n,1);for i=2:na(i-1,i-1)=pi^2/(12*n)+n;a(i-1,i)= pi^2/(24*n)-n;a(i,i-1)= pi^2/(24*n)-n;for j=1:nif j==i-1a(i,j)=a(i,i-1);else if j==ia(i-1,j-1)=2*a(i-1,i-1);else if j==i+1a(i,j)=a(i,i+1);elsea(i,j)=0;endendendendenda(n,n)=pi^2/(12*n)+n;for i=2:nf(i-1,i)=4/pi*cos((i-1)*pi/2/n)-8*n/(pi^2)*sin(i*pi/2/n)+8*n/(pi^2)*s in((i-1)*pi/2/n);endfor i=1:nf(i,i)=-4/pi*cos(i*pi/2/n)+8*n/(pi^2)*sin(i*pi/2/n)-8*n/(pi^2)*sin((i -1)*pi/2/n);end%b(j)=f(i-1,j)+f(i,j)for i=1:(n-1)b(i)=f(i,i)+f(i,i+1);endb(n)=f(n,n);tic;n=20;can=20;s=zeros(n^2,10);h=1/n;st=1/(2*n^2);A=zeros((n+1)^2,(n+1)^2);syms x y;for k=1:1:2*n^2s(k,1)=k;q=fix(k/(2*n));r=mod(k,(2*n));if (r~=0)r=r;else r=2*n;q=q-1;endif (r<=n)s(k,2)=q*(n+1)+r;s(k,3)=q*(n+1)+r+1;s(k,4)=(q+1)*(n+1)+r+1;s(k,5)=(r-1)*h;s(k,6)=q*h;s(k,7)=r*h;s(k,8)=q*h;s(k,9)=r*h;s(k,10)=(q+1)*h;elses(k,2)=q*(n+1)+r-n;s(k,3)=(q+1)*(n+1)+r-n+1;s(k,4)=(q+1)*(n+1)+r-n;s(k,5)=(r-n-1)*h;s(k,6)=q*h;s(k,7)=(r-n)*h;s(k,8)=(q+1)*h;s(k,9)=(r-n-1)*h;s(k,10)=(q+1)*h;endendd=zeros(3,3);B=zeros((n+1)^2,1);b=zeros(3,1);for k=1:1:2*n^2L(1)=(1/(2*st))*((s(k,7)*s(k,10)-s(k,9)*s(k,8))+(s(k,8)-s(k,10))*x+(s(k,9)-s(k,7))*y);L(2)=(1/(2*st))*((s(k,9)*s(k,6)-s(k,5)*s(k,10))+(s(k,10)-s(k,6))*x+(s (k,5)-s(k,9))*y);L(3)=(1/(2*st))*((s(k,5)*s(k,8)-s(k,7)*s(k,6))+(s(k,6)-s(k,8))*x+(s(k ,7)-s(k,5))*y);for i=1:1:3for j=i:3d(i,j)=int(int(((((diff(L(i),x))*(diff(L(j),x)))+((diff(L(i),y))*(dif f(L(j),y))))-((can^2)*L(i)*L(j))),x,0,1),y,0,1);d(j,i)=d(i,j);endendfor i=1:1:3for j=1:1:3A(s(k,(i+1)),s(k,(j+1)))=A(s(k,(i+1)),s(k,(j+1)))+d(i,j);endendfor i=1:1:3b(i)=int(int((L(i)),x,0,1),y,0,1);B(s(k,(i+1)),1)=B(s(k,(i+1)),1)+b(i);endendM=zeros((n+1)^2,n^2);j=n^2;for i=(n^2+n):-1:1if ((mod(i,(n+1)))~=1)M(:,j)=A(:,i);j=j-1;else continueendendpreanswer=M\B;answer=zeros((n+1)^2,1);j=1;for i=1:1:(n^2+n)if ((mod(i,(n+1)))~=1)answer(i)=preanswer(j);j=j+1;else answer(i)=0;endendZ=zeros((n+1),(n+1));for i=1:1:(n+1)^2s=fix(i/(n+1))+1;r=mod(i,(n+1));if(r==0)r=n+1;s=s-1;elseendZ(r,s)=answer(i);end[X,Y]=meshgrid(1:-h:0,0:h:1);surf(X,Y,Z);toc;t=toc;K=a ;B=b';U=inv(K)*Bfor i=1:nu(i,1)=4/(pi^2)*sin(pi*i/n/2);endue=U-u六、实验结果：程序中的变量can为问题中的k，为了方便比较，采用了画图的方式。