第8章 计算方法的MATLAB实现(MATLAB课件)

• yi=interp1(x,y,xi)返回在插值向量xi处的函 数向量yi，它是根据向量x和y插值而来。如 y是矩阵，则对y每一列进行插值，如xi中元 素不在x内，返回NaN。
• yi=interp1(y,xi)省略x，表示x=1:N，此时N 为向量y的长度或为矩阵y的行数。
2020/10/18
35
fx.m
• 首先编写上述非线性方程组的M文件fx.m
• function y=fx(x) • y(1)=0.1*(x(1)*x(1)+x(2)*x(2)+8); • y(2)=0.1*(x(1)*x(2)*x(2)+x(1)+8); • y=[y(1) y(2)];
ຫໍສະໝຸດ Baidu
2020/10/18
• x=
•1 •2 •3
2020/10/18
25
Gauss-Saidel迭代法
2020/10/18
26
gauss.m
• function s=gauss(a,b,x0,eps) • %gauss-seidel迭代法皆线性方程组 • %a为系数矩阵，b为方程组ax=b中的右边
的矩阵b，x0为迭代初值
20
Jacobi迭代法
2020/10/18
21
Jacobi.m
• function s=jacobi(a,b,x0,eps) • %jacobi迭代法皆线性方程组 • %a为系数矩阵，b为方程组ax=b中的右边
的矩阵b，x0为迭代初值 • if nargin==3 • eps=1.0e-6; • elseif nargin<3 • error • return • end
2020/10/18
28
程序实例
• 用上面编写的gauss函数求解下列方程组
10x1 2x2 x3 6 2x1 2x2 x3 10 x1 2x2 5x3 10
2020/10/18
29
程序实例
• >> a=[10 -2 -1;-2 2 -1;-1 -2 5]; • >> b=[6;10;10]; • >> x=gauss(a,b,[0 0 0]',eps)
• ans =
• 12 -3 3 • 0 3 -1 • 001
2020/10/18
12
程序实例
• >> triu(a,1) • ans = • 0 -3 3 • 0 0 -1 • 000 • >> triu(a,-1) • ans = • 12 -3 3 • -18 3 -1 • 011
2020/10/18
• x=
•4 • 13 •8
2020/10/18
30
8.3 非线性方程组数值解法
• 与线性方程组的求解一样，非线性方程组 的求解也是应用很广泛的课题。一般情况， 非线性方程组的数值解法主要采用迭代法 来求解。比较常用的迭代法主要有不动点 迭代法、Newton迭代法、拟Newton迭代法 等几种方法。
2020/10/18
33
• xx=fx(x);%第一次迭代 • while norm(xx-x)>=eps • x=xx; • xx=fx(x); • end • s=xx; • return
2020/10/18
34
程序实例
• 用不动点迭代法求解下面的方程组
x12 10x1x22 80 x1x22 x110x2 80
12x1 3x2 3x3 15 18x1 3x2 x3 15 x1 x2 x3 6
2020/10/18
9
程序实例
• >> a=[12 -3 3;-18 3 -1;1 1 1]; • >> b=[15;-15;6]; • >> x1=a\b • x1 = • 1.0000 • 2.0000 • 3.0000 • >> x2=inv(a)*b • x2 = •1 •2 •3
• ans = • -3 • -1 • >> diag(a,-1)
• ans = • -18 •1
2020/10/18
19
8.2.3 迭代解法
• 迭代解法非常适用于求解大型稀疏系数矩阵 的方程组，在线性方程组常用的迭代解法主 要有Jacobi迭代法、Gauss-Seidel迭代法。
2020/10/18
MATLAB 7.0从入 门到精通
哈尔滨工业大学(威海) 汽车工程学院
2020/10/18
1
课程主要内容
• 第1章 MATLAB简介 • 第2章 数值运算 • 第3章 单元数组和结构 • 第4章 字符串 • 第5章 符号运算 • 第6章 MATLAB绘图基础 • 第7章 程序设计 • 第8章 计算方法的MATLAB实现 • 第9章 优化设计 • 第10章 SIMULINK仿真初探
2020/10/18
2
第8章 计算方法的MATLAB实现
• 随着计算机的迅速发展与广泛运用，在众 多的领域，科学计算方法的应用越来越广 泛了，而MATLAB在进行科学计算方面有 着无与伦比的优势。本章介绍MATLAB在 计算方法中的运用。
2020/10/18
3
程序实例
• >> fzero('x^2-3*x+2',0) • ans = • 1.0000 • >> fzero('x^2-3*x+2',3) • ans = • 2.0000 • >> fzero('x^2-3*x+2',[0,4]) • ??? Error using ==> fzero • The function values at the interval
2020/10/18
10
8.2.2 线性方程组求解中的变换
• 上三角变换
• U=triu(x)返回矩阵x的上三角部分； • U=triu(x,k)返回第k条对角线以上部分的元
素。
2020/10/18
11
程序实例
• >> a=[12 -3 3;-18 3 -1;1 1 1]; • >> triu(a)
2020/10/18
40
程序实例
• 用上面编写的newton迭代函数求解下列 方程组
xx112x2210xx111x202x28800
2020/10/18
41
fx1.m和dfx1.m
• 首先，编写上述非线性方程组的M文件fx1.m • function y=fx1(x) • y(1)=x(1)*x(1)-10*x(1)+x(2)*x(2)+8; • y(2)=x(1)*x(2)*x(2)-10*x(2)+x(1)+8; • y=[y(1) y(2)]; • 然后，编写上述非线性方程组导数的M文件dfx1.m • function y=dfx1(x) • y(1)=2*x(1)-10; • y(2)=2*x(2); • y(3)=x(2)*x(2)+1; • y(4)=2*x(1)*x(2)-10; • y=[y(1) y(2);y(3) y(4)];
2020/10/18
7
8.2.1 直接解法
• 关于线性方程组的直接解法有许多种，比 如Gauss消去法、列主元消去法、平方根法 等。而在MATLAB中，线性方程组的直接 解法只需用符号“/”或“\”就解决问题。还 可以使用逆阵函数来求解： x=inv(A)*B。
2020/10/18
8
程序实例
• 求解下列方程组
初值，eps为允许误差
• if nargin==1 • eps=1.0e-6; • elseif nargin<1 • error • return • end
2020/10/18
39
• x1=fx1(x);%非线性方程组 • x2=-dfx1(x);%非线性方程组导数 • x3=inv(x2); • x0=x3*x1'; • while norm(x0)>=eps • x=x0'+x; • x1=fx1(x); • x2=-dfx1(x); • x3=inv(x2); • x0=x3*x1'; • end • s=x0'+x; • return
36
程序实例
• >> x=staticiterate([0 0]) • x= • 1.0000 1.0000
2020/10/18
37
Newton迭代法
2020/10/18
38
newtoniterate.m
• function s=newtoniterate(x,eps) • %newton迭代法解非线性方程组，x为迭代
2020/10/18
22
• D=diag(diag(a));%求出对角矩阵 • D=inv(D);%求出对角矩阵的逆 • L=tril(a,-1);%求出严格下三角矩阵 • U=triu(a,1);%求出严格上三角矩阵
• B=-D*(L+U); • f=D*b; • s=B*x0+f; • while norm(s-x0)>=eps • x0=s; • s=B*x0+f; • end • return
13
• 下三角变换
• U=tril(x)返回矩阵x的下三角部分； • U=tril(x,k)返回第k条对角线以上下部分的
元素。
2020/10/18
14
程序实例
• >> a=[12 -3 3;-18 3 -1;1 1 1]; • >> tril(a)
• ans =
• 12 0 0 • -18 3 0 • 111
2020/10/18
15
程序实例
• >> tril(a,1) • ans = • 12 -3 0 • -18 3 -1 • 111 • >> tril(a,-1) • ans = • 000 • -18 0 0 • 110
2020/10/18
16
• 对角变换
• U=diag(x)返回矩阵x主对角线上的元素，返 回结果是一列向量形式；
• if nargin==3 • eps=1.0e-6; • elseif nargin<3 • error • return • end
2020/10/18
27
• L=tril(a,-1);%求出严格下三角矩阵 • D=diag(diag(a));%求出对角矩阵 • U=triu(a,1);%求出严格上三角矩阵 • C=inv(D+L); • B=-C*U; • f=C*b; • s=B*x0+f; • while norm(s-x0)>=eps • x0=s; • s=B*x0+f; • end • return
• U=diag(x,k)返回第k条对角线上的元素值。 • 当x为向量时生成矩阵。
2020/10/18
17
程序实例
• >> a=[12 -3 3;-18 3 -1;1 1 1]; • >> diag(a)
• ans =
• 12 •3 •1
2020/10/18
18
程序实例
• >> a=[12 -3 3;-18 3 -1;1 1 1]; • >> diag(a,1)
2020/10/18
23
程序实例
• 用上面编写的jacobi函数求解下列方程组
10x1 2x2 x3 3 2x1 10x2 x3 15 x1 2x2 5x3 10
2020/10/18
24
程序实例
• >> a=[10 -2 -1;-2 10 -1;-1 -2 5]; • >> b=[3 15 10]'; • >> x=jacobi(a,b,[0 0 0]',eps)
2020/10/18
31
不动点迭代法
2020/10/18
32
staticiterate.m
• function s=staticiterate(x,eps) • %不动点迭代法解非线性方程组，x为迭代
初值，eps为允许误差 • if nargin==1 • eps=1.0e-6; • elseif nargin<1 • error • return • end
endpoints must differ in sign.
2020/10/18
6
8.2 线性方程组数值解法
• 线性方程组的求解不仅在工程技术领域涉 及到，而且在其他的许多领域也经常碰到， 因此这是一个应用相当广泛的课题。
• 关于线性方程组的数值解法一般分为两类： 一类是直接法，就是在没有舍入误差的情 况下，通过有限步四则运算求得方程组准 确解的方法。另一类是迭代法，就是先给 定一个解的初始值，然后按一定的法则逐 步求出解的各个更准确的近似值的方法。
2020/10/18
42
程序实例
• >> x=newtoniterate([0 0]) • x= • 11
2020/10/18
43
8.4 插值与拟合
• 在生产实践及科学实验中，插值与拟合的 应用非常广泛。下面，就对如何用 MATLAB来处理插值与拟合作一介绍。
2020/10/18
44
8.4.1 一维插值