第6章 MATLAB数值运算

x
求函数零点的函数。 x = fzero(fun,x0) fun: 函数的句柄、M文件函数名、匿名函数、 函数字符串。 x0:估计的零点位置。当x0为标量时，在它的 两侧寻找零点，当x0为一个二元向量[a,b]时， 在区间[a,b]之间寻找零点。 (1)估计的零点位置 将方程变换成h(x)=g(x)的形式，用ezplot 或fplot函数画出h(x)和g(x)的图象，二者的交 点即为解。
其解析解：
x
x n 1
1
3 xn h * f ( xn , t n )
e 30 t
x(1)=1/3;h=0.02; t=0:h:0.8; for i=1:length(t)-1 x(i+1)=x(i)+h*(-30*x(i)); end plot(t,x,'r*',t,x,'r-'),grid
fplot('-x/6',[-4*pi,4*pi]) grid hold off [m,n]=ginput
3 2
1
0
-1
-2
-3
-10
-5
0
5
10
function f=myfun(x) f=[sin(x(1))+x(2) x(1)+6*x(2)];
(2)求解
3
sin x y 0 x 6 y 0
i 1
ij
K
j
0.9
0.8
y c1 c2 e
t
0.7
0.6
下一步是用观测到的数据来估计 C=[c1 c2]。 建立方程：E*C=y
0 0.5 1 t 1.5 2 2.5
y
0.5
y c1 c2 e
E = 1.0000
1.0000
t
>>E = [ones(size(t)) exp(-t)]
1.0000
(3)四阶五级Runge-Kutta变步长算法(ode45 )
x n 1 x n K i hn f hn * i K i
i 1 6
t n i hn , x n i 1, 2 , ,6

j 1
xe 2 0
x
e
x
2/x 3
2 x
g ( x) 2 x
ezplot('exp(x)') hold on ezplot('2/x') hold off grid [x,y]=ginput(1)
2
1
h( x) e x
0
-1
-2
-3
-6
-4
-2
0 x
2
4
6
xe 2 0
x
(2)求解 函数字符串:
rank(A)=rand([A,B])=2
>>A=[1 –2 1 1;1 –2 3 –4;1 –2 1 5];B=[1;-1;5]; >>rank(A),rank([A,B]), X=A\B >> X=0 -0.0000 0 1.0000
6.1.2 解高次线性方程
x x 2x x 3 0
n
n 1
a n 1 x a n
将多项式的系数按降幂次序存放在一个行 向量中。缺少的幂次一定要补零。 p=[a0, a1,… an]
p( x) x 2 x 5
3
p=[1 0 -2 -5]
6.2.2 多项式乘除 多项式乘法实质上是系数向量的卷积 (convolution)，除法实质上是系数向量的解 卷积(deconvolution)。
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
n=3 12 10 8 6 4 2 0 -2 12 10 8 6 4 2 0 -2
n=5
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
6.2.4 多项式插值（Interpolation） 由已知数据点来估计其间的未知数据点。 一维插值函数： yi = interp1(x,y,xi,method) 二维插值函数： ZI = interp2(X,Y,Z,XI,YI,method)
a1=polyfit(x,y,1); x1=linspace(0,1,100); subplot(2,2,1) plot(x,y,'o',x1,polyval(a1,x1),'r') title('线性拟合')
n=1 12 10 8 6 4 2 0 -2 12 10 8 6 4 2 0 -2
n=2
y 0.476 0.3413 e t
0
0.5
1
1.5
2
2.5
(3)欠定方程(Underdetermined)，未知数个数多于方程个数， 当rank(A)=rand([A,B])，理论上有无穷个解，MATLAB可 求得一个基本解，其中至少有m个非零元素。
x1 2 x 2 x3 x 4 1 x1 2 x 2 x3 x 4 1 x 2x x 5x 5 2 3 4 1
a=[ 0.7800 0.5630 0.9130 0.6590] >>det(a) ans = 1.0000e-006 >>cond(a) ans = 2.1932e+006 >>s=svd(a) s= 1.48095205864320 0.00000067524130
(2)超定方程：实验数据的曲线拟合时，常遇到超定方程。 比如，实验中观测到下一组数据， t = [0 .3 .8 1.1 1.6 2.3]'; y = [.82 .72 .63 .60 .55 .50]'; t和y之间有什么关系呢？可用下面的函数近似描述：
>>x = fzero ( 'x*exp(x)-2' , [0 2] )

匿名函数:
>> fun=@ (x)x*exp(x)-2 >>x = fzero ( fun , [0 2] )

M文件函数名
function y=mfun(x) y=x*exp(x)-2 >> x = fzero ( mfun , [0 2] ) 得：x=0.8526
interp2
85
80
Temps
75
70
65
60 3 2 1 Depth 1 2 3 Width 4 5
6.3 常微分方程数值解法
微分方程数值解法是动态系统数字仿真的基础。
x f ( x, t ) x (t n 1 ) x (t n )

t nቤተ መጻሕፍቲ ባይዱ1
tn
f ( x , t ) dt
例：
width=1:5; depth=1:3; temps=[82 81 80 82 84; 79 63 61 65 81; 84 84 82 85 86]; di=1:0.2:3; wi=1:0.2:5; tc=interp2(width,depth,temps,wi,di','cubic'); mesh(wi,di,tc) hold on [ww,dd]=meshgrid(width,depth); plot3(ww,dd,temps,'ko') hold off
>>x=fsolve(@myfun,[2.8 -0.6]) x = 2.6788 -0.4465 >>x=fsolve(@myfun,[-2.8 0.6]) x =-2.6788 0.4465
2
1
0
-1
-2
-3
-10
-5
0
5
10
6 .2 多项式运算与插值法
6.2.1 多项式表示与创建
p ( x ) a0 x a1 x
0.7408
>>C=E\y C =0.4760 0.3413 C为观测数据的最小 二乘解。
1.0000
1.0000 1.0000 1.0000
0.4493
0.3329 0.2019 0.1003
y 0.476 0.3413 e
t
则：E*C=y，E的维数为6*2，为超定方程。
0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5
6.1.4 解非线性方程组
sin x y 0 x 6 y 0
函数： x = fsolve(fun,x0)，用来求解f(x)=0， 其中x是一个向量。 fun:关于向量x的M函数或匿名函数的句柄。 x0：寻找解的初始位置。 (1)估计解的位置
sin x y 0 fplot('-sin(x)',[-4*pi,4*pi]) x 6 y 0 hold on
5 4 3 2
>>p=[1 -1 2 1 0 3]; >>r=roots(p) r= 0.8679 + 1.1770i 0.8679 - 1.1770i -1.0000 0.1321 + 1.1770i 0.1321 - 1.1770i
6.1.3 解非线性方程
f ( x ) xe 2 0
把求微分变为求积分，然后用近似的方法求上 式中的积分，是数值积分法的基本思想。
2.3.1 欧拉法（Euler method）
用矩形面积近似曲边梯形的面积： f(x,t) 欧拉法递推公式如下：
x n 1 x n h * f ( x n , t n )
tn tn+1 t
误差
dx dt 30 x x (0) 1 / 3 0 t 0 .8
方程A*X=B,得X=A\B,如果A的维数是m*n则： m=n：恰定方程，可求得精确解 m>n：超定方程(Overdetermined)，可求得最小二乘解。 m<n：欠定方程(Underdetermined)，可求得一个基本解，
其中至少有m个非零元素。 (1)恰定方程 只要det(A)不等于0，即A非奇异时，方程有唯一的解， X=A\B或X=inv(A)*B，当矩阵的维数比较大时，前者的计算速 度和精度更高。 当A接近奇异时，要注意解的可靠性，可用函数cond(A) 来求A的条件数(Condition number)，条件数越接近1越好。
f(x,t)
误差
tn
tn+1
t
(2)四阶Runge-Kutta法 多取几个点，然后将其斜率加权平均得一等效斜 率，就得到四阶龙格库塔法。
x n 1 x n h * ( k1 2 k 2 2 k 3 k 4 ) / 6 k1 f ( x n , t n ) k 2 f ( x n k1 * h / 2, t n h / 2 ) k f ( x k * h / 2, t h / 2 ) n 2 n 3 k 4 f ( xn k 3 * h, t n h )
a(s) s 2s 3
2
b( s ) 4 s 5s 6
2
求a(s)*b(s): >>a = [1 2 3]; b = [4 5 6]; c = conv(a,b) c =4 13 28 27 18
6.2.3 多项式求值(Polynomial Evaluation)： >>p=[1 0 -2 -5] 3 p( x) x 2 x 5 >>polyval(p,5) ans = 110
6.3.2 龙格库塔算法(Runge-Kutta Method)
(1)二阶Runge-Kutta法 用梯形面积近似曲边梯形的面积：
x n 1 x n h * ( k1 k 2 ) / 2 k1 f ( x n , t n ) k f ( x k * h, t h ) n 1 n 2
6.2.3 多项式拟合（Fitting） p = polyfit(x,y,n) x,y是观测到的数据点向量，长度为N，n为多项式 的阶次，p为多项式的系数向量，长度为n+1。拟合 的算法是最小二乘法。 例：x=0 : 0.1 : 1;
y=[-0.47,1.97,3.28,6.16,7.08,7.34,7.66,8.36,9.28,9.60,10.2];
MATLAB语言
电气学院 叶满园 myye@4y.com.cn
第6章 MATLAB数值运算

解析解:符号运算工具箱(Symbolic Math Toolbox) 数值解:很多工程问题不存在或很难得到解析解、解析解
计算太复杂时，需要用数值解法。
6.1 解方程 6.1.1 解线性方程组 方程A*X=B的解为：X=inv(A)*B=A\B， A\B称为A左除B，左除时要求两矩阵行数相等。 方程X*A=B的解为：X=B*inv(A)=B/A， A/B称为A右除B，右除时要求两矩阵列数相等