Matlab关于数值计算的实现要点

Matlab关于数值计算的实现

摘要：数值计算（numerical computation computation），主要研究更好的利用计算机更好的进行数值计算，解决各种数学问题。数值分析包括离散傅里叶变换，考虑截断误差，计算误差，函数的敛散性与稳定性等。在数学方面，数值计算的主要研究数值微分与积分，数据的处理与多项式计算，最优化问题，线性方程与非线性方程的求解，常微分方程的数值求解等。同时，数值计算在物理，化学，经济等方面也有研究，本文暂且不表。M atlab软件历经二十多年来的发展，已成为风靡世界的数学三大软件（matlb，Mathematica l，Maple）之一，在数学类科技应用软件中在数值计算方面首屈一指。Matlab以矩阵为数据操作的基本单位，使得矩阵运算十分便捷快速，同时Matlab还提供了海量的计算函数，而且使用可靠地算法进行计算，能使用户在繁复的数学运算中解脱，Matlab还具有方便且完善的图形处理功能，方便绘制二维和三维图形并修饰。

目录

1.数值分析(离散傅里叶变换，考虑截断误差，计算误差，函数

的敛散性与稳定性)

2.数值计算(数值微分与积分，数据的处理与多项式计算，

最优化问题，线性方程与非线性方程的求解，常微分方程的数值求解)

3.图形处理功能(方便绘制二维和三维图形并修饰)

4.总结

1.数据统计与分析

Matlab 可以进行求矩阵的最大最小元素，平均值与中值，关于矩阵元素的求和与求积，累加和与累乘积，标准方程，相关系数，元素排序。现在以求标准方差举例说明Matlab 的实现。

在Matlab 中，实现标准方差计算的函数为std 。对于向量(Y ），std （Y ）实现返回一个标准方差，而对于矩阵（A ），std （A ）返回一个行向量，该行向量的每个元素对应着矩阵A 各行或各列的标准方差。一般调用std 函数的格式为std （A ，flag ，dim ）

Dim 取1或者2分别对应求各列或各行的标准方差，flag 取1时，按照标准方差的计算公式

∑-=-=N

i x x S i N 1

2

1)(11来计算。若flag 取2，则用公式

∑-==N

i x x S i N

1

2

2)

(1

进行计算。默认的flag 取值为0，dim 取值为1。课本page143

2. 离散傅里叶变换

离散傅里叶变换广泛应用于信号的分析，光谱和声谱分析、全息技术等各个领域。但直接计算dft 的运算量与变化的长度N 的平方成正比，当N 较大时，计算量太大。随着计算机技术的迅速发展，在计算机上进行离散傅里叶变换计算成为可能。特别是快速傅里叶变换算法的出现，为傅里叶变换的应创造了条件。

（1）：傅里叶变换算法的简述。

傅立叶变换是一种分析信号的方法，它可分析信号的成分，也可用这些成分合成信号。许多波形可作为信号的成分，比如正弦波、方波、锯齿波等，傅立叶变换用正弦波作为信号的成分. f(t)是t 的周期函数，如果t 满足狄里赫莱条件：在一个以2T

为周期内f(X)连续或只有有限个第一类间断点，附f （x ）单调或可划分成有限个单调区间，则F （x ）以2T 为周期的傅里叶级数收敛，和函数S （x ）也是以2T 为周期的周期函数，且在这些间断点上，函数是有限值；在一个周期内具有有限个极值点；绝对可积。则有下图①式成立。称为积分运算f(t)的傅立叶变换，

②式的积分运算叫做F(ω)的傅立叶逆变换。F(ω)叫做f(t)的像函数，f(t)叫做

F(ω)的像原函数。F(ω)是f(t)的像。f(t)是F(ω)原像。 ①傅立叶变换

②傅立叶逆变换

(2)离散傅里叶变化的实现

Matlab提供了对向量或对矩阵进行离散傅里叶变换的函数.

1. g(x):返回向量的离散傅里叶变换.设晓得长度为n，若为2的幂次，则为以2为基数快速

傅里叶变换，否则为运算速度很慢的非2次幂次的算法。对于矩阵x，个g（x）应用于矩阵每一列。

2.g（x，n）计算n点离散傅里叶变化。限定向量的长度为n，若x的长度小于n，则不足

部分不上0；若大于n，则删去超出的那些元素。对于矩阵X，它同样应用于矩阵的每一列，只是限定了向量长度为n

3.g（x，[],dim）:这是对于矩阵而言的函数调用格式，，前者的功能与g（x）

基本相同，而后者则与g（x，n）基本相同。只是当参数dim=1时，改函数作用于x的每一列；当dim=2时，则作用x的每一行。

值得一提的是，当一只给出的样本数n0是2的幂次时，可以取一个n使他大于n0且是2的幂次，然后利用函数格式g（x，n）进行快速的傅里叶变换。是七速度大大的提升。

例如：

给定函数：

g（t）=12sin(2pi*10t+pi\4)+5cos(2pi*40t)

取N=100，试从x为0-1采样，绘制相应的振幅-频率图。

程序如下：

N=100 采样点数

T=1; % 采样时间终点

t=linspace(0,T,N);

% 给出N个采样时间

ti(I=1:N)

x=12*sin(2*pi*10*t+pi/4)+5*cos(2*pi*40*t);

% 求各采样点样本值x

dt=t(2)-t(1);

% 采样周期f=1/dt;

% 采样频率(Hz)

X=fft(x); % 计算x的快速傅立叶变换X

F=X(1:N/2+1); % F(k)=X(k)(k=1:N/2+1)

f=f*(0:N/2)/N; % 使频率轴f从零开始

plot(f,abs(F),'-*') % 绘制振幅-频率图

xlabel('Frequency');

ylabel('|F(k)|');

2.数值计算

一、数值微分与积分

Matlab数值积分

MATLAB提供最简单的积分函数是梯形法trapz，我们先说明梯形法语法trapz(x,y)，其中x,y 分别代表数目相同的阵列或矩阵，而y与x的关系可以由是一函数型态（如y=sin(x)）或是不以函数描述的离散型态。

我们看一简单积分式

以下为MA TLAB 的程式

（and(1,101)产生101个均匀随机数,每个数都介于0-1之间）

>> x=0:pi/100:pi;

>> y=sin(x);

>> k=trapz(x,y)

k =

1.9998

clear

x=sort(rand(1,101)*pi);

y=sin(x);

trapz(x,y)

ans =

1.9981

2.TLAB 另外提供二种积分函数，它们分别是辛普森法quad 和牛顿-康兹法quad8。三种方法的精确度由低而高，分别为trapz, quad, quad8。

由于这二种方法依据的积分法不同于梯形法，因此它们的语法就和trapz 不同；其语法为quad('function',a,b) （quad8语法相同），其中function是一已定义函数的名称（如sin, cos, sqrt, log 等），而a, b是积分的下限和上限。和trapz比较，quad, quad8不同之处在于这二者类似解析式的积分式，只须设定上下限及定义要积分的函数；而trapz则是针对离散点型态的数据做积分。

我们看一简单积分式

以下为MATLAB 的程式

>> a=0; b=0.5;

>> kq=quad('sqrt',a,b)

kq =

0.2357

>> kq8=quad8('sqrt',a,b)

kq8 = 0.2357

?--2

03521

dx x x

clear

fun=inline('1./(x.^3-2*x-5)') ezplot(fun,[0,2])

[q,n]=quad(fun,0,2) % [q,n] = quad(fun,a,b,?) %同时返回

%函数计算的次数n q =

-0.4605 n = 53

3.积分与三重积分

??+ππ

π

2)]cos()sin([dxdy

y x x y

fun=inline('y*sin(x)+x*cos(y)') Q=dblquad(fun,pi,2*pi,0,pi) Q = -9.8698

[x,y]=meshgrid(pi:.1:2*pi,0:.1:pi); z=fun(x,y); mesh(x,y,z)

Matlab数值微分

若有一函数f(x) 在x=a 的微分可表示为

上述微分项在几何上的意义为在x=a 点的切线斜率. 下面介绍的差分表示式，数值差分即是用来求数值微分的方法。

差分表示法

我们知道一微分项的计算，可以在二相邻点x+h 和x 间函数取下列极限求得

若将原连续的空间以多个离散点取代。则上述的极限以离散点的方式计算，即是以下的差分式(difference equation)

上式被称为前向差分

而高阶微分项可以利用低阶微分项来计算，举例来说一个二阶微分式可以表示为

上述提及的后向差分式，在MA TLAB 有对应的diff 函数来计算二相邻点的差值，它的语法为diff(x)，其中x代表一组离散点。假设有x, y(x)的数据为

x=[1 3 5 7 9], y=[1 4 9 16 25]

则diff(x)=[2 2 2 2], diff(y)=[3 5 7 9]，

注意二者皆以后向差分计算且数据点只剩4 个而不是5个。而dy/dx 的数值微分则为

dy/dx=diff(y)./diff(x)。

因此要计算下列多项式在[-4, 5] 区间的微分

>> x=linspace(-4,5); % 产生100个x的离散点

>> p=[1 -3 -11 27 10 -24]; %被积函数各项的系数

>> f=polyval(p,x);

>> plot(x,f) % 将多项式函数绘图

>> title('Fifth-deg. equation')

>> dfb=diff(f)./diff(x); % 注意要分别计算diff(f)和diff(x) >> xd=x(2:length(x)); % 注意只有99个df值，而且是对应x2,x3,...,x100的点

>> plot(xd,dfb ) % 将多项式的微分项绘图

>> title('Derivative of fifth-deg. equation')

>> num=f(3:length(f))-f(1:length(f)-2); % 注意中央差分是f(k+1)-f(k-1)

>> deno=x(3:length(f))-x(1:length(f)-2); % 注意中央差分是x(k+1)-x(k-1)

>> df_c=num./deno;

>> xd=x(2:length(x)-1); % xd的点数只有98个

>> plot(xd,df_c)

>> title('Derivative of fifth-deg. polynomial')

二、线性方程与非线性方程的求解

1.线性方程组求解

在MATLAB中，线性方程的求解分为两种：一是直接解法，就是在没有舍入误差的情况下，通过有限的矩阵的初等运算来求的方程组的解；二是符号解法。就是先顶一个解的初值，然后按照一定的秩代法运算进行逐步的不仅。求出更精确的近似解。

（1）：直接解法

线性方程的直接解法多数基于高斯消元法、主元素消元法、平方根法和追赶法。只需要调用相应的函数或运算符即可完成线性方程的求解。最简单的方法那就是左除法。如

用直接法求线性方程组的解法

X1+2X2+3X3+4X4=10

X1+3X2+4X3+5X4=9

X1-3X2-4X3+5X4=6

X1-5X2+3X3-4X4=11

命令如下：

A=[1,2,3,4;1,3,4,5;1,-3,-4,5;1,-5,3,-4]

B=[10,9,6,11]

X=A\B

(2)矩阵分解求线性方程

A=[

5,1,-1,1;

1,3,-1,-1;

-1,-3,-1,5;

1,1,-1,4];

b=[1;2;3;4];

[l,u]=lu(A); %LU分解

X=u\(L\b)

（3）矩阵求逆法

A=[

5,1,-1,1;

1,3,-1,-1;

-1,-3,-1,5;

1,1,-1,4];

b=[1;2;3;4];

X=inv（A）*b

2.非线性方程组求解

方程求根的一般形式是求下列方程的根：f ( x ) = 0 实际上，就是寻找使函数f ( x ）等于零的变量x ，所以求方程（l ）的根， 也叫求函数

f ( x ）的零点。如果变量x 是列阵，则方程就代表方程组。当方程中的函数f (x ）是有限个指数、对数、三角、反三角或幂函数的组合时，则方程（ l ）被称为超越方程，例如 e-x - sin （π x/ 3) ＋lnx = 0 就是超越方程

一、符号方程求解

在MATLAB 中，求解用符号表达式表示的代数方程可由函数solve 实现，其调用格式为：

solve(s)：求解符号表达式s 的代数方程，求解变量为默认变量。 当方程右端为0时，方程可以不标出等号和0，仅标出方程的左端。 solve(s,v)：求解符号表达式s 的代数方程，求解变量为v 。

solve(s1,s2,…,sn,v1,v2,…,vn)：求解符号表达式s1,s2,…,sn 组成的代数方程组，求解变量分别v1,v2,…,vn 。

例1． 解下列方程。

1．2

2

144212-+=-++x x x x x= solve('1/(x+2)+4*x/(x^2-4)=1+2/(x-2)', 'x')

2．17433=---x x x f=sym('x-(x^3-4*x-7)^(1/3)=1') x= solve(f)

3．143sin 2=??? ?

?

-πx

x= solve('2*sin(3*x-pi/4)=1') 4．010=-+x xe x

x= solve('x+x*exp(x)-10', 'x') %仅标出方程的左端

2 迭代法

迭代法是计算数学中的一种重要方法，用途很广，求解线性方程组和矩阵特征值时也要用到它。这里结合非线性方程的迭代法求解，介绍一下它的基本原理。

迭代法基本原理

迭代法的基本原理就是构造一个迭代公式，反复用它得出一个逐次逼近方程根的数列，数列中每个元素都是方程根的近似值，只是精度不同。

迭代法求解方程

f ( x ) = 0 （1） 时，先把方程等价地变换成形式

f ( x ) = x －g(x) = 0 ， （2） 移项得出：

x = g(x) ( 3 )

若函数g (x ）连续，则称（3）为迭代函数。用它构造出迭代公式：

x k+1= g ( x k ） , k = 0 , l , 2 ， … ( 4 ) 从初始值 x 0出发，便可得出迭代序列：

{ x k ｝= x 0, x 1, x 2,….x k ,….. （ 5 ） 如果迭代序列（5 ）收敛，且收敛于x*，则由式（4）有：

()()()()()0***lim 1==-=-+∞

→x f x x g x x g k k

可见 x*便是方程（l ）的根。

3.图形处理功能

Matlab 提供了一系列的绘图函数，用户不需要过多的考虑绘图的细节，只需要给出一些基本参数就能得到所需图形，这类函数称为高层绘图函数。此外，Matlab 还提供了直接对图形句柄进行操作的低层绘图操作。这类操作将图形的每个图形元素（如坐标轴、曲线、文字等）看做一个独立的对象，系统给每个对象分配一个句柄，可以通过句柄对该图形元素进行操作，而不影响其他部分。绘图可分为两种，是二维图形绘图和三维图形绘图。

一、 二维图形绘图

Matlab 最基本的绘图函数为plot 函数，，用于绘制xy 平面上的线性坐标曲线图，提供新坐标和对应的y 坐标，即plot （x ，y ）。

如 在[0 , 2pi]区间，绘制曲线y=2*exp(-0.5*x).*sin(2*pi*x);

程序如下：在命令窗口中输入以下命令 >> x=0:pi/100:2*pi;

>> y=3*exp(-0.7*x).*sin(4*pi*x); >> plot(x,y)

2． 含多个输入参数的plot 函数

plot 函数可以包含若干组向量对，每一组可以绘制出一条曲线。含多个输入参数的plot 函数调用格式为：plot(x1，y1，x2，y2，…，xn ，yn) 如下列命令可以在同一坐标中画出3条曲线。 >> x=linspace(0,2*pi,50);

>> plot(x,sin(x),x,5*sin(x),x,7*sin(x))

3．含选项的plot函数

Matlab提供了一些绘图选项，用于确定所绘曲线颜色。这些选项如表所示：

线型颜色

- 实线b蓝色

: 虚线g绿色

-. 点划线r红色

-- 双划线c青色

m品红

y黄色

k黑色

w白色

例用不同颜色在同一坐标内绘制曲线及其包络线。

>> x=(0:pi/100:2*pi)';

>> y1=3*exp(-0.5*x)*[1,-1];

>> y2=3*exp(-0.5*x).*sin(2*pi*x);

>> x1=(0:12)\3;

>> y3=3*exp(-0.5*x1).*sin(2*pi*x1);

>> plot(x,y1,'k:',x,y2,'b--',x1,y3,'g');

分别绘制了黑、蓝、绿三种曲线。

二．绘制三维曲线的基本函数

最基本的三维图形函数为plot3，它将二维绘图函数plot的有关功能扩展到三维空间，可以用来绘制三维曲线。其调用格式为：

plot3（x1，y1，z1，选项1，x2，y2，z2，选项2，…）

其中每一组x，y，z组成一组曲线的坐标参数，选项的定义和plot的选项一样。当x，y，z 是同维向量时，则x，y，z对应元素构成一条三维曲线。当x，y，z是同维矩阵时，则以x，y，z对应列元素绘制三维曲线，曲线条数等于矩阵的列数

如绘制空间曲线

该曲线对应的参数方程为

t=0:pi/100:2*pi;

x=5*cos(t);

y=2*sqrt(2)*sin(t);

z=-2*sqrt(2)*sin(t);

plot3(x,y,z,'p');

title('Line in 3-D Space');

text(0,0,0,'origin');

xlabel('X');ylabel('Y');zlabel('Z');grid

三维绘图不做详细的介绍。

4.总结

总的来说，matlab对于我们来说，提供了很大的方便，它的功能很强大，除具备卓越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。对于数学专业来说。我们不必执着于用手来计算一切，我们只需通过一些简单的带代码，即可将繁杂的的数学积分、微分、方程组的求解，数据的处理等很快的计算出来。同事还能够绘制二三维图形，其立体型很具体，生动形象。因此，对于我们数学专业的人来说，学好matlab这一门课程是很必要的。