matlab程序设计实践.doc

matlab程序设计实践

———————————————————————————————— 作者：
———————————————————————————————— 日期：

17

个人收集整理 勿做商业用途
MATLAB程序设计实践
1、编程实现以下科学计算算法,并举一例应用之。(参考书籍《精
通ＭＡＬＡＢ科学计算》,王正林等著，电子工业出版社，２００９
年)
“里查森迭代法线性方程组求解"
解：
算法说明：
里查森迭代法是最简单的迭代法，它的迭代公式为：xk+1=(I-A）*xk+b;在MATLAB中编程实现的里查森迭代法函数为：richason。
功能:用里查森迭代法求线性方程组
调用格式：[x，n]=richason（A,b，x0,eps，M）
其中,A为线性方程组的系数矩阵；
b为线性方程组的常数向量；
x0为迭代初始向量;
eps为解的精度控制(此参数可选）；
M为迭代步数控制(此参数可选）；
x为线性方程组的解;
n为求出所需精度的解实际的迭代步数.
里查森迭代法的MATLAB程序代码如下：
function [x,n］ = richason(A,b,x0，eps,M)
％采用里查森迭代法求线性方程组Ax=b的解
％线性方程组的系数矩阵：A
％线性方程组的常数向量：b
％迭代初始向量：x0
%解的精度控制：eps
%迭代步数控制：M
%线性方程组的解:x
％求出所需精度的解实际的迭代步数：n
if(nargin==3）
eps=1。0e-6； %eps表示迭代精度
M=200； %M表示迭代步数的限制值
elseif（nargin==4)
M=200；
end
I=eye(size(A））;
x1=x0;
x=(I—A）＊x0+b;
n=1；
%迭代过程
while(norm（x—x1）〉eps）
x1=x;
x=（I-A）*x1+b;
n=n+1； ％n为最终求出解时的迭代步数
if(n〉=M）
disp（’Warning:迭代次数太多,可能不收敛！’);
return;
end
end

实例:用里查森迭代法求以下线性方程组,其中初始值取为[0 0 0]

输入：
〉〉 A=[1.0170 -0.0092 0。0095；
—0.0092 0。9903 0.0136；
0.0095 0。0136 0。9898］;
>> b=［1 0 1］';
>> x0=［0 0 0］’;
〉> ［x，n]=richason(A，b,x0)
输出的计算结果为:
x =

0。9739
-0。0047
1。0010

输出的迭代次数为:
n =

5

经过5步迭代，理查森迭代法求出了方程的解为：
[x1，x2，x3］=［0.9738,—0。0047,1。0010］

对上述迭代计算结果进行验证，在MATLAB命令窗口中输入如下程序：
〉〉 A＊x
输出结果为：
ans =
1。0000
0。0000
1。0000
经检验，计算结果正确。
程序运算截图如下：

开始
流程图:
否
、源程序
n=1
是
Warning:迭代次数太多，可能不收敛
结束

否
是
x1=x;
n=n+1
否
输出结果
读取数据
nargin==3？
eps=10--6 ，最大步数M=200
最大步数M为200
x=(I-A)*x0+b
norm(x-x1)>eps?
是
n>=200?
例题流程图
输入系数矩阵A

↓
输入初始向量x0
及常数向量b

[x,n]=richason(A,b,x0)

↓

输出计算结果

↓

输出迭代次数

↓

A*x验证结果

↓

解：
（1）算法说明
分析已给方程可知,为拉普拉斯方程，在MATLAB工具箱PDETOOL中可看成椭圆型方程，转化为标准形式如下：

因此，对应的c=—1,a=0,f=0，然后根据给出的边界约束条件，在微分方程工具箱中选择所需要的条件,
① Dirichlet条件

② Neumann条件

其中n是上的单位外法矢量，g，q，h和r是定义在上的函数。
（题目中Γ1与Γ2分别代表x+y=2与x—y=2这两条边界线）

（2)操作流程
设置坐标限
选择Options栏中Axes Limits选项，输入坐标范围

绘制区域图
点击绘制多边形键 画出要求的区域图

设置边界条件
选择Boundary中的Boundary Mode，设置为边界模式双击各条边界线，由方程组中已知边界条件设定
设置方程参数
点击 ，将已知方程对照标准偏微分方程形式,知c=-1，a=0，f=0。

剖分网格
按顺序点击 两按钮,细分网格.

绘制温度分布图

点击绘制三维示意图:

(3)简易流程图
开始
绘制要求区域图
设置边界条件
设置方程参数
剖分网格
绘制温度分布示意图
结束

实验1 用GUI方式解下列PDE

解：
（1）算法说明
同上题，由已给方程可知,为拉普拉斯方程，在PDETOOL中可看成椭圆型方程，转化为标准形式如下：

因此,对应的c=1,a=0,f=0,然后根据给出的边界约束条件,在微分方程工具箱中选择所需要的条件，
③ Dirichlet条件

④ Neumann条件

其中n是上的单位外法矢量，g，q,h和r是定义在上的函数。

（2）操作流程

设置坐标限

绘制区域图

设置边界条件
u｜x=0=y（3-y）

u｜y=0=sin（π/4*x)

设置方程参数

划分网格

绘制特征值对应的函数图形
二维图形：

三维图形:

（3）简易流程图

开始
绘制要求区域图
设置边界条件
设置方程参数
剖分网格
绘制特征值对应的函数图形
结束