第7章计算方法的MATLAB实现讲稿

第7章 计算方法的MATLAB实现

第7章 计算方法的MATLAB实现

计算方法主要研究数学问题的数值解，涉及的内容广。利用MATLAB提供的部分函数解决，可以实现某些情况下的数值求解。本章作简要介绍。

7.1 一元非线性方程求解

求解一元非线性方程的方法主要有二分法、割线法、牛顿法等。本节主要介绍两个可以求解一元非线性方程的函数，即fzero函数和roots函数。 7.1.1 fzero函数

用fzero函数求一元非线性方程的零点。该函数的简单调用格式为： ●x=fzero(fun,x0)：如果x0为标量，则试图寻找函数fun在x0附近的零点。fun参数是一个M文件函数或匿名函数的函数句柄。函数fzero返回的值x靠近fun函数改变符号的位置。如果搜索失败，则返回NaN。搜索区间扩展到发现inf、NaN或复数值时终止搜索。如果x0是一个长度为2的向量，则fzero函数假设x0是一个区间，其中fun(x0(1))与fun(x0(2))异号。如果符号相同，则会出错。给出相应函数值异号的区间可以保证fzero函数返回一个fun函数改变符号的位置附近的值。

●[x,fval]=fzero(…)：返回解x和解x处目标函数fun的值。

●[x,fval,exitflag]=fzero(…)：还返回一个exitflag值，描述fzero函数的退出条件。返回值及其描述如下所示：

>0，表示函数找到了零值点x； ?0，表示没有发现零值点。

例7－1 求方程(x?20)(x?1)(x?5)(x?30)?0的解。 程序代码： clear;clc;

f=@(x)(x+20)*(x+1)*(x-5)*(x-30);%创建匿名函数 x1=fzero(f,-50)%解一元非线性方程 x2=fzero(f,1) x3=fzero(f,6) x4=fzero(f,50) 运行结果： x1 =

- 133 -

第7章 计算方法的MATLAB实现

-20 x2 = -1 x3 = 5 x4 = 30

7.1.2 roots函数

用roots函数计算多项式的根，该函数的调用格式为： ●r=roots(c)：返回一个列向量，其元素为多项式c的解。 例7－2 求多项式方程x?5x?x?5?0的解。 程序代码：

clear;clc;

P=[1 -5 1 -5];%表示多项式 r=roots(P)%求多项式方程的解

P1=poly(r)%由多项式的根返回多项式的系数 运行结果：

r =

5.0000 -0.0000 + 1.0000i -0.0000 - 1.0000i P1 =

1.0000 -5.0000 1.0000 -5.0000

注：poly和roots互为逆函数，函数poly返回多项式的系数。

327.2 非线性方程组的数值解法

用MATLAB求解线性方程组的解，在6.3中已介绍了基于矩阵变换的直接解法，除此方法之外，还有很多方法，比如，Jocabi迭代法、Gauss-Seidel迭代法和SOR(超松驰)迭代法等，这里不作详细介绍。可参阅有关文献。实际工程中得到的数学模型往往具有非线性的特点，得到其解析解比较困难。一般采用迭代法求解非线性方程组。比较常见的迭代方法有不动点迭代法、Newton迭代法和拟Newton迭代法等几种。本节仅介绍不动点迭代法。

设含有n个未知数和n个方程的非线性方程组记为：

- 134 -

第7章 计算方法的MATLAB实现

F(x)?0

其中，x为由n个未知数构成的向量，F为由n个函数构成的向量。

将方程组F(x)?0改写为便于迭代的等价形式：x??(x)，由此得出不

动点迭代法的迭代公式：xk?1??(xk)。

k*k??k*如果得到的序列{x}满足limx?x，则x就是?的不动点。这样就可

以求出线性方程组的解（近似解）。

据此，编写不动点迭代法的M文件(文件名为example7_2)： function s=example7_2(x,eps)

%用不动点迭代法求非线性方程组的解 %x为迭代初值，eps为允许误差值。 if nargin==1 eps=1e-6; elseif nargin<1 error return end

x1=example7_2a(x);% example7_2a是函数x??(x)的文件名。 while norm(x1-x)>=eps x=x1;

x1=example7_2a(x);%循环迭代 end s=x1; return

2?x12?11x1?x2?9?0例7－3 用不动点迭代法求方程组?的解。 2xx?x?6x?8?012?122?x1??(x12?x2?9)/11解：变形为? 2x?(xx?x?8)/6121?2首先创建函数代码（文件名为example7_2a.m）： function y=example7_2a(x)

y(1)=(x(1)*x(1)+x(2)*x(2)+9)/(-11); y(2)=(x(1)*x(2)*x(2)+x(1)+8)/6; 再作图形，程序代码： clear;clc;

- 135 -

第7章 计算方法的MATLAB实现

ezplot('x1^2+11*x1+x2^2+9',[-12,6]); hold on

ezplot('x1*x2^2+x1-6*x2+8',[-12,6]); hold off

title('x1^2+11*x1+x2^2+9和x1*x2^2+x1-6*x2+8的图形')%修改标题 运行结果见图7－1：

图7－1

最后调用函数example7_2求方程组解： x=example7_2 ([0,0]) 运行结果：

x=

-1.0000 1.0000

2?x12?11x1?x2?9?0 这样就可以求得方程组?的一个解，从图7－1可2?x1x2?x1?6x2?8?02?x12?11x1?x2?9?0以看出，方程组?还有一个解，为了求得它，还需要对2?x1x2?x1?6x2?8?0原方程组作另外变形，这里不再作进一步的讨论。

- 136 -

第7章 计算方法的MATLAB实现

7.3 插值

插值计算在数据拟合和数据平滑等方面应用普遍。MATLAB提供了用最近邻插值、线性插值、三次样条插值、三次插值和FFT插值法进行一维、二维、三维和高维插值的函数。 7.3.1 一维插值

MATLAB中有两种一维插值，即多项式插值和基于FFT的插值。 ⒈ 多项式插值

函数interp1进行一维插值。一维插值是进行数据分析和曲线拟合的重要手段。interp1函数使用多项式技术，用多项式函数拟合所提供的数据，并计算目标插值点上的插值函数值，其最常用的语法形式是：

yi=interp1(x,y,xi,method)

x和y为给定数据的向量，长度相同。xi为包含要插值的点的构成的向量，method是一个可选的字符串，指定一种插值方法，包括：

●最近邻插值（method='nearest'）：该方法将插值点的值设置为已知数据点中距离最近的点的值；

●线性插值（method='linear'）：该方法用线性函数拟合每对数据点，并返回xi处的相关函数值；

●三次样条插值（method='spline'）：该方法用三次样条函数拟合每对数据点，用spline函数在插值处进行三次样条插值；

●三次插值（method='pchip'或'cubic'）：该方法用pchip函数对向量x和y进行分段三次Hermite插值。

这几种方法在速度、内存和平滑性等方面有所差别，使用时可以根据需要进行选择，包括：

◣最近邻插值是最快的方法，但是，利用它得到的结果平滑性最差； ◣线性插值比最近邻插值要占用更多的内存，运行时间略长。与最近邻法不同，它生成的结果是连续的，但是在顶点处有坡度变化；

◣三次样条插值的运行时间相对来说最长，内存消耗比三次插值略少。它生成的结果平滑性最好。但是，如果输入数据不很均匀，可能会得到意想不到结果；

◣三次插值需要更多内存，而且运行时间比最近邻插值和线性插值的长。但是，使用此法时，插值数据及其导数都是连续的。

例7－4 程序代码：

- 137 -