《偏微分方程数值解法》课程设计

题    目：     中心差分解两点边值问题               

*****************************************

学    院：          理学院                    

专    业：       信息与计算科学               

班    级：          0708032                     

学    号：         *********                 

指导老师：***

2010年12月25日

一．题目

用中心差分格式计算如下两点边值问题

已知其精确解为

二．理论

作为模型，考虑两点边值问题：

…………（1.1）

…………（1.2）

假定是给定的常数。

1.建立差分格式

（1）.区域网格剖分

首先取个节点：

将区间分成个小区间：

于是得到区间的一个网格剖分。记，称为网格最大步长。用表示网格内点，，，的集合，表示内点和界点的集合。

取相邻节点的中点，称为半整数点。则由节点

又构成的一个网格剖分，称为对偶剖分。

（2）.微分方程的离散，建立相应差分格式

用差商代替微商，将方程（1.1）在内点离散化.注意对充分光滑的，由Taylor展式有

………（1.3）

 ………（1.5）

由（1.5）减（1.4），并除以，得

…………（1.6）

令则由（1.3）（1.6）知，边值问题的解满足方程：

 …………（1.7）

其中

 …………（1.8）

为差分算子的截断误差，舍去，便得逼近边值问题（1.1）（1.2）的差分方程：

 …………（1.9）

i=1,2,…,N-1，

由方程（1.7）（1.9），截断误差可表示为

 …………（1.10）

当网格均匀，即时差分方程（1.9）简化为

 …………（1.11）

这相当于用一阶中心差商，二阶中心差商依次代替（1.1）的一阶微商和二阶微商的结果。这个方程就是中心差分格式。

截断误差为：

 …………（1.12）

所以截断误差按或的阶为。

在本题中，，，， 

，

因为r=0方程（1.11）的系数对角矩阵是三对角矩阵。我们可以用消元法或迭代法求解方程组（1.1）（1.2）

式（1.11）用方程组展开：

写成矩阵形式为：

2.收敛性分析

根据（1.10）我们引进误差

则误差函数满足下列差分方程：

于是收敛性及收敛速度的估计问题，就归结到通过右端（截断误差）估计误差函数的问题。

由（1.12）我们知，有

从而差分方程满足相容条件。

若引进记号

，，，

，，

设

则可将（1.9）改写为

将差分解表成

 …………（2.1）

其中满足

 …………（2.2）

而满足

 …………（2.3）

先估计，由

 …………（2.4）

据差分格林公式

再利用柯西不等式，有常数使

 …………（2.5）

将不等式（2.6）用于（2.5）右端，则

 …………（2.6）

解差分方程（2.2，易得）

从而

这样，

 …………（2.7）

利用范数，从（2.7）推出

 …………（2.8）

因为

因此

 …………（2.9）

联结（2.1）（2.7）及（2.9）即得差分解的先验估计：

 …………（2.10）

其中

不等式（2.10）说明差分解连续依赖于右端和边值，因此差分格式（1.11）关于右端及边值稳定.

根据定理1.1 ：  若边值问题的解u充分光滑，差分方程按满足相容条件且关于右端稳定，则差分解按收敛到边值问题的解，且有和相同的收敛阶。所以差分方程的解的收敛速度为。

三．程序代码：

clc

clf

clf

syms x;

a区间界点

b区间界点

p这是p函数

q这是q函数

f=-exp(x)*(2*x+1)+(sin(x)+1+x)*x*(x-1);%这是f函数

r这是r函数.

N=10将区间划分的等分,这里控制！

h这里确定步长

value_of_f=zeros(N-1,1);%这是f

diag_0=zeros(N-1,1);%确定A的对角元

diag_1=zeros(N-2,1);%确定A的偏离对角的上对角元

diag_2=zeros(N-2,1);%确定A的偏离对角的下对角元

X=a:h:b;

u边界条件

u边界条件

for j=2:N

e获取对角元素

for j=3:N

e获取A的第三条对角

for j=2:N-1

e获取A的第二条对角

for j=2:N;

e获取F值

value_of_f(1)=value_of_f(1)+u_a*(subs(p,{x},{(X(2)+X(1))/2}))/(h^2);

value_of_f(N-1)=value_of_f(N-1)+u_b*(subs(p,{x},{(X(N)+X(N+1))/2}))/(h^2);

A=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1);%组装系数矩阵

format long

U差分解

%fprintf('%11.5f',U)

fprintf('\\n');

dx=X(2:N);

p精确解

%fprintf('%11.5f',precise_value)

d误差

deta_max=max(abs(deta));%最大误差

fprintf('最大的误差是%f\\n',deta_max)

plot(X(2:N),U,'b*',X(2:N),precise_value,'r--')  %差分解与精确解对比表

figure();

p误差图

结果：

   X的值

N=20,h=0.05(图中蓝点和紫线分别是精确值图和差分图)

四．结论：

由本题可以总结出，通过用中心差分解两点边值所得的数值解能够较好地逼近方程的精确解，且区域剖分得越细，即步长越小，数值解与精确解的误差就越小，数值解越逼近精确解。下载本文