实验09 数值微积分与方程数值解（第6章）

《数学软件》课内实验

王平

实验09 数值微积分与方程数值求解

（第6章 MATLAB数值计算）

一、实验目的

1. 掌握求数值导数和数值积分的方法。 2. 掌握代数方程数值求解的方法。 3. 掌握常微分方程数值求解的方法。 二、实验内容

1. 求函数在指定点的数值导数

xf(x)?1程序及运行结果： x2x36x2x3x2,x?1,2,3

022. 用数值方法求定积分

(1) I1??2?0cost2?4sin(2t)2?1dt的近似值。

程序及运行结果：

1

(2) I2? 2??0ln(1?x)dx

1?x2程序及运行结果：

3. 分别用3种不同的数值方法解线性方程组

?6x?5y?2z?5u??4?9x?y?4z?u?13? ??3x?4y?2z?2u?1??3x?9y?2u?11程序及运行结果： 4. 求非齐次线性方程组的通解

?2x1?7x2?3x3?x4?6??3x1?5x2?2x3?2x4?4 ?9x?4x?x?7x?2234?1程序及运行结果（提示：要用教材中的函数程序line_solution）： 5. 求代数方程的数值解

(1) 3x+sinx-ex=0在x0=1.5附近的根。

程序及运行结果（提示：要用教材中的函数程序line_solution）：

(2) 在给定的初值x0=1，y0=1，z0=1下，求方程组的数值解。

?sinx?y2?lnz?7?0?y3 ?3x?2?z?1?0?x?y?z?5?0?程序及运行结果： 2

6. 求函数在指定区间的极值

x3?cosx?xlogx(1) f(x)?在(0,1)内的最小值。

ex程序及运行结果： 332(2) f(x1,x2)?2x1在[0,0]附近的最小值点和最小值。 ?4x1x2?10x1x2?x2程序及运行结果： 7. 求微分方程的数值解，并绘制解的曲线

?xd2ydy?5?y?0?dx2dx?? ?y(0)?0?y'(0)?0???

程序及运行结果（注意：参数中不能取0，用足够小的正数代替）：

令y2=y,y1=y '，将二阶方程转化为一阶方程组：

1?'5y?y??1x1xy2?' ?y2?y1?y(0)?0,y(0)?02?1? 8. 求微分方程组的数值解，并绘制解的曲线

?y'1?y2y3?y'??yy?213 ?y'??0.51yy12?3??y1(0)?0,y2(0)?1,y3(0)?1

程序及运行结果： 3

三、实验提示

四、教程：第6章 MATLAB数值计算（2/2）

6.2 数值微积分 p155 6.2.1 数值微分

1. 数值差分与差商

对任意函数f(x)，假设h>0。

? 向前差分：?f(x)?f(x?h)?f(x) ? 向后差分：?f(x)?f(x)?f(x?h)

? 中心差分：?f(x)?f(x?h/2)?f(x?h/2) 当步长h充分小时，有

?f(x) h?f(x)? 向后差商：f'(x)?

h?f(x)? 中心差商：f'(x)?

h? 向前差商：f'(x)?2. 数值微分的实现

MATLAB没有提供求数值导数的函数，只有计算向前差分的函数diff。 ? DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,…,n-1。 ? DX=diff(X,n)：计算X的n阶向前差分。 例如，diff(X,2)=diff(diff(X))。

? DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列；dim=2，按行。

例6.18 （向前差分）求1~3阶差分p156

设x由[0,2π]间均匀分布的10个点组成，求sinx的1~3阶差分。

4

X=linspace(0,2*pi,10); Y=sin(X); DY=diff(Y) D2Y=diff(Y,2) D3Y=diff(Y,3) DY = 0.6428 0.3420 -0.1188 -0.5240 -0.6840 -0.5240 -0.1188 0.3420 0.6428 D2Y = -0.3008 -0.4608 -0.4052 -0.1600 0.1600 0.4052 0.4608 0.3008 D3Y = -0.1600 0.0556 0.2452 0.3201 0.2452 0.0556 -0.1600 例6.19 （数值导数）3种方法求导p157

f(x)?x3?2x2?x?12?6x?5?5x?2

用不同的方法求f(x)的数值导数，并在同一个坐标系中做出f '(x)的图像。 f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2'); g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2-x+12)/2+1/6./(x+5).^(5/6)+5');%导函数 x=-3:0.01:3; p=polyfit(x,f(x),5);%用5次多项式p拟合f(x) dp=polyder(p); %对拟合多项式p求导数dp dpx=polyval(dp,x);%求dp在假设点的函数值 dx=diff(f([x,3.01]))/0.01;%直接对f(x)求数值导数 gx=g(x);%求函数f的导函数g在假设点的导数 plot(x,dpx,x,dx,'.',x,gx,'-'); 5

6.2.2 数值积分 p157

1. 数值积分基本原理 求解定积分的数值方法： ? 梯形法

? 辛普生(Simpson)?法

? 牛顿－柯特斯(Newton-Cotes)法 基本思想：

将整个积分区间[a,b]分成n个子区间 [x, x]，i=1,2,…,n，其中x=a，xn+1=b

i

i+1

1

这样求定积分问题就分解为求和问题。 2. 数值积分的实现

(1) 被积函数是一个解析式

调用格式： quad(fname,a,b,tol,trace)

quadl(fname,a,b,tol,trace)

? fname是被积函数名。

? a和b分别是定积分的下限和上限。 -6

? tol用来控制积分精度，默认tol=10。

? trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认trace=0。

6

例6.20（解析式函数数值积分）两种方法p158

用两种不同方法求 I??e01?x2dx

I = 0.746824180726425 I = 0.746824133988447 J = 0.746824180726425 J = 0.746824133988447 format long I=quad('ex',0,1) I=quadl('ex',0,1) g=inline('exp(-x.^2)');%内联函数 J=quad(g,0,1) %无 ' J=quadl(g,0,1) format short function ex=ex(x) ex=exp(-x.^2); (2) 被积函数由一个表格定义 调用格式：trapz(X,Y) ? Y=f(X)

? X=[x1,x2,...,xn], xi

例6.21（表格函数数值积分）p159

用trapz函数计算 I?format long; X=0:0.01:1; Y=exp(-X.^2); trapz(X,Y) format short; (3) 二重积分数值求解

?10e?x2dx

ans = 0.746818001467970 I??dc?baf(x,y)dxdy

调用格式：

I=dblquad(f,xmin,xmax,ymin,ymax,tol,trace) 例6.22（二重数值积分）p160

I??1?1?2?2e?x/2sin(x2?y)dxdy

2format long global ki; ki=0; I=dblquad('fxy',-2,2,-1,1) ki f=inline('exp(-x.^2/2).*sin(x.^2+y) ', 'x', 'y'); J=dblquad(f,-2,2,-1,1) format short 7

function f=fxy(x,y) global ki; ki=ki+1;%调用次数 f=exp(-x.^2/2).*sin(x.^2+y); I = 1.574493189744944 ki = 1050 J = 1.574493189744944 表 数值微积分函数和命令 p155~160

函数/命令 diff inline quad trapz dblquad 定义内联函数 解析式函数数值积分 表格定义的函数数值积分 二重积分数值求解 说 明 向前差分的函数（用于数值微分） 6.3 离散傅立叶变换 p160 6.3.1 离散傅立叶变换算法简述

在某时间片等距抽取N个抽样时间tm处的样本值f(tm)，记f(m)，m=1,2,...,N。 ? 称向量F(k)(k=1,2,...,N)为f(m)的一个离散傅里叶变换，公式：

F(k)??f(m)e?j2?(m?1)(k?1)/N,k?1,2,m?1N,N

又称F(k)为f(m)的离散频谱。

? 由F(k)逆求f(m)的过程，称离散傅立叶逆变换，公式：

f(m)??F(k)ej2?(m?1)(k?1)/N,m?1,2,k?1N,N

6.3.2 离散傅立叶变换的实现 p161

一维离散傅立叶变换函数：

(1) fft(X)：返回向量X的离散傅立叶变换。 ? 设X的元素个数为N。

? 若N为2的幂次，则为以2为基数的快速傅立叶变换。 ? 否则为运算速度很慢的非2幂次的算法。 ? 对于矩阵X，它应用于矩阵的每一列。 (2) fft(X,N)：计算N点离散傅立叶变换。 ? 它限定向量的长度为N。

8

? 若X的长度小于N，则不足部分补零。 ? 若大于N，则删去超出N的那些元素。 ? 对于矩阵X，它应用于矩阵的每一列。

(3) fft(X,[ ],dim)或fft(X,N,dim)：前者与fft(X)基本相同，后者与fft(X,N)基本相同。 ? dim=1时，作用于X的每一列。 ? dim=2时，作用于X的每一行。

注：当已知给出的样本数N0不是2的幂次时，可取一个N使它大于N0且是2的幂次，然后利用函数格式fft(X,N)或fft(X,N,dim)便可进行快速傅立叶变换。

一维离散傅立叶逆变换函数，调用格式：

(1) ifft(F)：返回F的一维离散傅立叶逆变换。 (2) ifft(F,N)：为N点逆变换。

(3) ifft(F,[ ],dim)或ifft(F,N,dim)：则由N或dim确定逆变换的点数或操作方向。

例6.23 （快速傅里叶变换）p161

给定数学函数

x(t)=12sin(2π×10t+π/4)+5cos(2π×40t)

取N=128，试对t从0~1秒采样，用fft作快速傅立叶变换，绘制相应的振幅-频率图。

在0~1秒时间范围内采样128点，确定采样周期和采样频率。由于离散傅立叶变换时的下标应是从0到N-1，故在实际应用时下标应该前移1。

对离散傅立叶变换，其振幅|F(k)|是关于N/2对称的，故只须使k从0到N/2即可。 N=128;%采样点数 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)|') 9

ix=real(ifft(X));%求逆变换，结果只取实部 plot(t,x,t,ix,':')%逆变换结果和原函数的曲线 norm(x-ix)%逆变换结果和原函数之间的距离 ans = 2.2403e-014 10

表 离散傅立叶变换函数和命令 p160~163

函数/命令 fft ifft 说 明 一维离散傅立叶变换函数 一维离散傅立叶逆变换函数 6.4 线性方程组求解 p163 6.4.1 直接解法

1. 利用左除运算符的直接解法

对于线性方程组Ax=b，可利用左除运算符“\\”求解： x=A\\b

? A为方阵，若是奇异的，则给出警告信息。 ? b为列向量，得一个解。 ? b为矩阵，得多个解。

例6.24（左除）求解线性方程组p163

?2x1?x2?5x3?x4?13?x?5x?7x??9?124 ?2x?x?x?6?234??x1?6x2?x3?4x4?0 11

A=[2,1,-5,1; 1,-5,0,7; 0,2,1,-1; 1,6,-1,-4]; b=[13,-9,6,0]'; x1=A\\b x2=inv(A)*b x1 = -66.5556 25.6667 -18.7778 26.5556 x2 = -66.5556 25.6667 -18.7778 26.5556 2. 利用矩阵的分解求解线性方程组

矩阵分解指按一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。 常见的矩阵分解： ? LU分解 ? QR分解

? Cholesky分解 ? Schur分解

? Hessenberg分解 ? 奇异分解 (1) LU分解

将一个矩阵表示为一个下三角阵和一个上三角阵的乘积形式。 已经证明，只要方阵A是非奇异的，LU分解是可行的。 调用格式：

? [L,U]=lu(A)：产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，满足A=LU。

? [L,U,P]=lu(A)：产生一个上三角阵U和一个下三角阵L及一个置换矩阵P，满足PA=LU。

线性方程组Ax=b的解：

x=U\\(L\\b) 或 x=U\\(L\\P*b)

注意：使用第一种格式时，矩阵L往往不是一个下三角矩阵，但可通过行交换成一个下三角阵。

例 LU分解 p164

?1?11??

A??5?43???11??2?A=[1,-1,1; 5,-4,3; 2,1,1]; [L,U]=lu(A)%(1) LU=L*U [L,U,P]=lu(A)%(2) L = 0.2000 -0.0769 1.0000 1.0000 0 0 0.4000 1.0000 0 U = 5.0000 -4.0000 3.0000 12

LU=L*U inv(P)*L*U 0 2.6000 -0.2000 0 0 0.3846 LU = 1 -1 1 5 -4 3 2 1 1 L = 1.0000 0 0 0.4000 1.0000 0 0.2000 -0.0769 1.0000 U = 5.0000 -4.0000 3.0000 0 2.6000 -0.2000 0 0 0.3846 P = 0 1 0 0 0 1 1 0 0 LU = 5 -4 3 2 1 1 1 -1 1 ans = 1 -1 1 5 -4 3 2 1 1 例6.25（LU分解）求解线性方程组p166

?2x1?x2?5x3?x4?13?x?5x?7x??9?124 ??2x2?x3?x4?6??x1?6x2?x3?4x4?0A=[2,1,-5,1; 1,-5,0,7; 0,2,1,-1; 1,6,-1,-4]; b=[13,-9,6,0]'; [L,U]=lu(A); x1=U\\(L\\b) [L,U ,P]=lu(A); x2=U\\(L\\P*b) (2) QR分解

x1 = -66.5556 25.6667 -18.7778 26.5556 x2 = -66.5556 25.6667 -18.7778 26.5556 13

A为方阵，调用格式：

? [Q,R]=qr(A)：产生一个正交矩阵Q和一个上三角矩阵R，满足A=QR。

? [Q,R,E]=qr(A)：产生一个正交矩阵Q、一个上三角矩阵R及一个置换矩阵E，满足AE=QR。

线性方程组Ax=b的解：

x=R\\(Q\\b) 或 x=E*(R\\(Q\\b))

例 QR分解 p166

?1?11??

A??5?43????2710??>> A=[1,-1,1; 5,-4,3; 2,7,10]; >> [Q,R]=qr(A) %(1) Q = -0.1826 -0.0956 -0.9785 -0.9129 -0.3532 0.2048 -0.3651 0.9307 -0.0228 R = -5.4772 1.2780 -6.5727 0 8.0229 8.1517 0 0 -0.5917 >> QR=Q*R QR = 1.0000 -1.0000 1.0000 5.0000 -4.0000 3.0000 2.0000 7.0000 10.0000 >> [Q,R,E]=qr(A) %(2) Q = -0.0953 -0.2514 -0.9632 -0.2860 -0.9199 0.2684 -0.9535 0.3011 0.0158 R = -10.4881 -5.4347 -3.4325 0 6.0385 -4.2485 0 0 0.4105 E = 0 0 1 0 1 0 1 0 0 >> Q*R/E ans = 1.0000 -1.0000 1.0000 14

5.0000 -4.0000 3.0000 2.0000 7.0000 10.0000 >> 例6.26（QR分解）求解线性方程组p167

?2x1?x2?5x3?x4?13?x?5x?7x??9?124 ?2x?x?x?6?234??x1?6x2?x3?4x4?0 (3) Cholesky分解 A是对称正定矩阵，调用格式：

? R=chol(A)：产生一个上三角阵R，使R'R=A。若A为非对称正定，则输出出错信息。

? [R,p]=chol(A)：不输出出错信息。当A为对称正定的，则p=0，R与上述格式得到的结果相同；否则p为一个正整数。若A为满秩，则R为一个阶数为q=p-1的上三角阵，且满足R'R=A(1:q,1:q)。

线性方程组Ax=b的解：

x=R\\(R'\\b)

例Cholesky分解p168（对称正定）

11??2?

A??12?1????1?13?? 15

例6.27（Cholesky分解）求解线性方程组p168（非对称正定）

?2x1?x2?5x3?x4?13?x?5x?7x??9?124 ??2x2?x3?x4?6??x1?6x2?x3?4x4?0 16

6.4.2 迭代解法 p168

迭代解法适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括： ? Jacobi迭代法

? Gauss-Serdel迭代法 ? 超松弛迭代法 ? 两步迭代法

1. Jacobi迭代法

线性方程组Ax=b，若A为非奇异方阵，即aii≠0(i=1,2,…,n)，则可将A分解为 A=D-L-U

? D为对角阵，其元素为A的对角元素。 ? L与U为A的下三角阵和上三角阵。 于是Ax=b化为：

x=D-1(L+U)x+D-1b

与之对应的Jacobi迭代公式为：

x(k+1)=D-1(L+U)x(k)+D-1b

若序列{x(k+1)}收敛于x，则x必是Ax=b的解。

Jacobi迭代法的函数文件Jacobi.m p170

function [y,n]=jacobi(A,b,x0,eps) if nargin==3 eps=1.0e-6; elseif nargin<3 error return end D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 B=D\\(L+U); f=D\\b; y=B*x0+f; n=1; %迭代次数 while norm(y-x0)>=eps x0=y; y=B*x0+f; n=n+1; end 例6.28（Jacobi迭代法）求解线性方程组 p170

设迭代初值为0，迭代精度为10-6。

17

?10x1?x2?9???x1?10x2?2x3?7 ??2x?10x?623?>> A=[10,-1,0;-1,10,-2;0,-2,10]; >> b=[9,7,6]'; >> [x,n]=jacobi(A,b,[0,0,0]',1.0e-6) x = 0.9958 0.9579 0.7916 n = 11 2. Gauss-Serdel迭代法

迭代公式：

x(k+1)=(D-L)-1Ux(k)+(D-L)-1b

和Jacobi迭代相比，Gauss-Serdel迭代用新分量代替旧分量，精度高。

Gauss-Serdel迭代法的函数文件gauseidel.m p171

function [y,n]=gauseidel(A,b,x0,eps) if nargin==3 eps=1.0e-6; elseif nargin<3 error return end D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 G=(D-L)\\U; f=(D-L)\\b; y=G*x0+f; n=1; %迭代次数 while norm(y-x0)>=eps x0=y; y=G*x0+f; n=n+1; end 例6.29（Gauss-Serdel迭代法）求解线性方程组 p171

设迭代初值为0，迭代精度为10-6。

18

?10x1?x2?9???x1?10x2?2x3?7 ??2x?10x?623?>> A=[10,-1,0;-1,10,-2;0,-2,10]; >> b=[9,7,6]'; >> [x,n]= gauseidel (A,b,[0,0,0]',1.0e-6) x = 0.9958 0.9579 0.7916 n = 7 例6.30（两种迭代法）求解线性方程组 p171

设迭代初值为0。是否收敛？

?12?2??x1??9??111??x???7? ???2????1??22???x3????6??>> a=[1,2,-2;1,1,1;2,2,1]; >> b=[9;7;6]; >> [x,n]=jacobi(a,b,[0;0;0]) x = -27 26 8 n = 4 >> [x,n]=gauseidel(a,b,[0;0;0]) x = NaN NaN NaN n = 1012 6.4.3 求线性方程组的通解 p172

Ax=b，n为未知变量个数。 ① rank(A)=n，有唯一解： x=A\\b ② b=0，为齐次方程组，x=0为平凡解。

rank(A)

19

③ b≠0

? rank(A)=rank([A,b])=n，有唯一解

x=A\\b 或 x=pinv(A)*b

? rank(A)=rank([A,b])

通解=特解+齐次方程组的基础解系

用A\\b求特解。

? rank(A)

求线性方程组的函数文件line_solution.m p172

function [x,y]=line_solution(A,b) [m,n]=size(A); y=[ ]; if norm(b)>0 %非齐次方程组 if rank(A)==rank([A,b]) if rank(A)==n disp('有唯一解x'); x=A\\b; else disp('有无穷个解，特解x，基础解系y'); x=A\\b; y=null(A,'r'); end else disp('无解'); x=[ ]; end else %齐次方程组 disp('有零解x'); x=zeros(n,1); if rank(A)

?x1?2x2?3x3?x4?1??3x1?x2?5x3?3x4?2 ?2x?x?2x?2x?334?12>> A=[1,-2,3,-1; 3,-1,5,-3; 2,1,2,-2]; >> b=[1;2;3]; >> [x,y]=line_solution(A,b) 无解 20

x = [ ] y = [ ] 例6.32（求线性方程组的通解）p174

?x1?x2?3x3?x4?1??3x1?x2?3x3?4x4?4 ?x?5x?9x?8x?0234?1>> format rat >> A=[1,1,-3,-1; 3,-1,-3,4; 1,5,-9,-8]; >> b=[1;4;0]; >> [x,y]=line_solution(A,b) 有无穷个解，特解x，基础解系y Warning: Rank deficient, rank = 2, tol = 8.8373e-015. > In line_solution at 11 x = 0 0 -8/15 3/5 y = 3/2 -3/4 3/2 7/4 1 0 0 1 表 线性方程组求解函数和命令 p163~174

函数/命令 x=A\\b x=inv(A)*b lu qr chol Jacobi迭代公式 Gauss-Serdel迭代公式 null 说 明 左除求解线性方程组 对方阵LU分解 对方阵QR分解 对称正定矩阵的Cholesky分解 x(k+1)=D-1(L+U)x(k)+D-1b x(k+1)=(D-L)-1Ux(k)+(D-L)-1b 求齐次方程组的基础解系 21

6.5 非线性方程与最优化问题求解 p174 6.5.1 非线性方程数值求解

1. 单变量非线性方程求解 调用格式：

z=fzero(filename,x0,tol,trace)

? filename是待求根的函数文件名。

? x0为搜索的起点。函数可能有多个根，但fzero只给出离x0最近的那个根。 ? tol控制结果的相对精度，缺省tol=eps。

? trace?迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省trace=0。

例6.33（单变量非线性方程求解）p175

f(x)?x?>> fzero('fz',-5) ans = -5.1926 >> fzero('fz',1) ans = 0.1926 function y=fz(x) y=x-1/x+5; 2. 非线性方程组F(X)=0的求解 调用格式：

x=fsolve(filename,x0,option)

? x为返回的解。

? filename是用于定义需求解的非线性方程组的函数文件名。 ? x0是求根过程的初值。

? option为最优化工具箱的选项设定。

最优化工具箱提供了20多个选项，用户可使用optimset命令将它们显示出来。 常用选项： ① Display：决定函数调用时中间结果的显示方式，'off'不显示，'iter'每步都显示，'final'

只显示最终结果。 ② LargeScale：是否用大规模问题算法，取值'off'或'on'。 ③ Maxlter：最大允许迭代次数，选择空矩阵取默认值400。 ④ TolFun：目标函数误差容限，选择空矩阵取默认值10-6。 ⑤ TolX：自变量误差容限，选择空矩阵取默认值10-6。 可用option=optimset命令调入一组默认选项值。

若想改变其中某个选项，则可调用optimset()函数完成：

optimset('Display', 'off');

1?5 x在x0=-5和x0=1作为迭代初值时的零点。 22

或用结构体属性的方式设置新参数：

option=optimset;

option.LargeScale='off';

例6.34 （非线性方程组求解） p176

求非线性方程组在(1,1,1) 附近的解并对结果进行验证。

?sinx?y?z2ex?0? ?x?y?z?0?xyz?0?function F=myfun(X) x=X(1); y=X(2); z=X(3); F(1)=sin(x)+y+z^2*exp(x); F(2)=x+y+z; F(3)=x*y*z; >> X=fsolve('myfun',[1,1,1]',optimset('Display','off')) X = 0.0224 -0.0224 -0.0000 >> q=myfun(X) q = 1.0e-006 * -0.5931 -0.0000 0.0006 例6.35 （非线性方程组求解） p176

求圆和直线的两个交点。 圆： x2?y2?z2?9

?3x?5y?6z?0 直线： ?x?3y?6z?1?0?该问题即为求解方程组：

?x2?y2?z2?9?0??3x?5y?6z?0 ?x?3y?6z?1?0?分析：

使用fsolve函数求解方程组时，必须先估计出方程组的根的大致范围。 所给直线的方向数是(-12,24,-14)，故其与球心在坐标原点的球面的交点大致是(-1,1,-1)和(1,-1,1)。以这两点作为迭代初值。 function F=fxyz(X) 23

x=X(1); y=X(2); z=X(3); F(1)=x^2+y^2+z^2-9; F(2)=3*x+5*y+6*z; F(3)=x-3*y-6*z-1; >> X1=fsolve('fxyz',[-1,1,-1],optimset('Display','off')) X1 = -0.9508 2.4016 -1.5259 >> X2=fsolve('fxyz',[1,-1,1],optimset('Display','off')) X2 = 1.4180 -2.3361 1.2377 6.5.2 无约束最优化问题求解 p177

无约束最优化问题的一般描述为：

minf(x)

x其中x=[x1, x2, ..., xn]'，该数学表示的含义亦即求取一组x，使得目标函数f(x)为最小。

故这样的问题又称为最小化问题。 调用格式：

① [x,fval]=fminbnd(fname,x1,x2,options)

求一元函数在(xl,x2)区间中的极小值点x和最小值fval。 ② [x,fval]=fminsearch(fname,x0,options)

基于单纯形算法求多元函数的极小值点x和最小值fval。 ③ [x,fval]=fminunc(fname,x0,options)

基于拟牛顿法求多元函数的极小值点x和最小值fval。

例6.36 （无约束最优化问题求解）一元函数 p177

求函数：

f(x)?x?1?5 x在区间(-10,1)和(1,10)上的最小值点。 function f=fx(x) f=x-1/x+5; >> [x1,fmin1]=fminbnd('fx',-10,-1) x1 = -9.9999 fmin1 = -4.8999 >> f=inline('x-1/x+5'); >> [x2,fmin2]=fminbnd(f,-10,-1) x2 = -9.9999 24

fmin2 = -4.8999 >> fminbnd(f,1,10) ans = 1.0001 例6.37 （无约束最优化问题求解）多元函数 p178

设

y2z22f(x,y,z)?x???

4xyz求的f在(0.5,0.5,0.5)附近的最小值。 function f=fxyz(u) x=u(1); y=u(2); z=u(3); f=x+y.^2./x/4+z.^2./y+2./z; >> [U,fmin]=fminsearch('fxyz',[0.5,0.5,0.5]) U = 0.5000 1.0000 1.0000 fmin = 4.0000 6.5.3 有约束最优化问题求解 p178

xs.t.G(x)?0minf(x)

其中x=[x1, x2, ..., xn]'，求一组x，使目标函数f(x)最小，且满足约束条件G(x)≤0。

约束条件：

① 线性不等式约束：Ax≤b ② 线性等式约束：Aeqx=beq ③ 非线性不等式约束：Cx≤0 ④ 非线性等式约束：Ceqx=0 ⑤ x的下界和上界：Lbnd≤x≤Ubnd 调用格式：

[x,fval]=fmincon(fname,x0,A,b,Aeq,beq,Lbnd,Ubnd,NonF,options) ? x、fval、fname、x0和options的含义与求最小值函数相同。 ? 其余参数为约束条件。

? NonF为非线性约束函数的M文件名。 ? 若某个约束不存在，用空矩阵表示。

25

例6.38 （有约束最优化问题求解）p179

min2f(x)?0.4x2?x12?x2?x1x2?13x130?x1?0.5x2?0.4?s..t?0.5x1?x2?0.5?x?0,x?02?1

function f=fop(x) f=0.4*x(2)+x(1)^2+x(2)^2-x(1)*x(2)+1/30*x(1)^3; x0=[0.5;0.5]; A=[-1,-0.5;-0.5,-1]; b=[-0.4;-0.5]; lb=[0;0]; option=optimset; option.LargeScale='off'; option.Display='off'; [x,f]=fmincon('fop',x0,A,b,[ ],[ ],lb,[ ],[ ],option) Warning: Options LargeScale = 'off' and Algorithm = 'trust-region-reflective' conflict. Ignoring Algorithm and running active-set method. To run trust-region-reflective, set LargeScale = 'on'. To run active-set without this warning, use Algorithm = 'active-set'. > In fmincon at 412 In Untitled2 at 8 x = 0.3400 0.3300 f = 0.2456 表 非线性方程与最优化问题求解函数p174~179

函数 fzero fsolve optimset fminbnd fminsearch fminunc fmincon 说 明 单变量非线性方程求解 非线性方程组的求解 选项设置/显示命令 求一元函数在区间中的极小值 单纯形算法求多元函数的极小值 拟牛顿法求多元函数的极小值 有约束下求极小值 26

6.6 常微分方程的数值求解 p179

常微分方程的初值问题：

y'=f(t,y), t0≤t≤T y(t0)=y0

数值解法：求y(t)在节点t0

y0,y1,...,ym称为常微分方程初值问题的数值解。 一般取ti=t0+ih, i=0,1,...,m，h称为步长。 常用的数值解法： ? 欧拉(Euler)法

? 龙格—库塔(Runge-Kutta)法 ? 线性多步法 ? 预报校正法

6.6.1 龙格－库塔法简介

龙格—库塔公式：

hy(t0?ih)?yi?yi?1?(k1?2k2?2k3?k4)6其中，

k1?f(ti?1,yi?1)hhk2?f(ti?1?,yi?1?k1)22

hhk3?f(ti?1?,yi?1?k2)22k4?f(ti?1?h,yi?1?hk3)6.6.2 龙格－库塔法的实现

调用格式：

[t,y]=ode23(fname,tspan,y0) [t,y]=ode45(fname,tspan,y0)

? fname是定义f(t,y)的函数文件名，函数须返回一个列向量。 ? tspan形式为[t0,tf]，表示求解区间。 ? y0是初始状态列向量。

? t和y分别给出时间向量和相应的状态向量。

第一种格式采用二阶、三阶龙格—库塔法。 第二种格式采用四阶、五阶龙格—库塔法。

例6.39 （常微分方程的数值求解）一阶p181

设有初值问题：

27

?y2?t?2,0?t?1?y'? 4(t?1)??y(0)?2?试求其数值解，并与精确解相比较（精确解为y(t)?t?1?1）。 function y=funt(t,y) y=(y^2-t-2)/4/(t+1); t0=0;tf=10; y0=2; [t,y]=ode23('funt',[t0,tf],y0);%求数值解 y1=sqrt(t+1)+1; %求精确解 plot(t,y,'b.',t,y1,'r-')

例6.40 （常微分方程的数值求解）二阶p181

已知一个二阶线性系统的微分方程为：

?d2x?2?ax?0,a?0 ?dt?x(0)?0,x'(0)?1?其中a=2，绘制系统的时间响应曲线和相平面图。

令x2=x,x1=x '，二阶方程转化为一阶方程：

'?x2?x1?' ?x1??ax2?x(0)?0,x(0)?11?2function xdot=sys(t,x) 28

xdot=[-2*x(2);x(1)]; t0=0;tf=20; [t,x]=ode45('sys',[t0,tf],[1,0]); [t,x] subplot(1,2,1); plot(t,x(:,2));%解曲线，即t-x subplot(1,2,2); plot(x(:,2),x(:,1));%相平面曲线，即x-x' axis equal ans = 0 1.0000 0 0.0001 1.0000 0.0001 0.0001 1.0000 0.0001 0.0002 1.0000 0.0002 … 19.7511 -0.9422 0.2352 19.8340 -0.9747 0.1556 19.9170 -0.9937 0.0738 20.0000 -0.9991 -0.0090

表 常微分方程的数值求解函数p179~183

函 数 ode23 ode45 说 明 二阶、三阶龙格—库塔法 四阶、五阶龙格—库塔法 29

6.7 稀疏矩阵 p183 6.7.1 矩阵存储方式

MATLAB的矩阵有两种存储方式： 1．完全存储方式

将矩阵的全部元素按列存储。 2．稀疏存储方式

仅存储矩阵所有的非零元素的值及其位置，即行号和列号。 稀疏存储方式也是按列存储的。

例（矩阵存储方式）p183

?1000??

A??0500????2007??完全存储方式（按列）： 1，0，2，0，5，0，0，0，0，0，0，7 稀疏存储方式（按列）： (1,1),1，(3,1),2，(2,2),5，(3,4),7 6.7.2 稀疏存储方式的产生

1. 完全存储方式转化为稀疏存储方式

A=sparse(S)

将矩阵S转化为稀疏存储方式的矩阵A。

当矩阵S是稀疏存储方式时，则函数调用相当于A=S。

例6.41（将完全存储方式转化为稀疏存储方式）p183

设

?2?0?A??0??0??00000000051000000?0??0? ??1??5??A = (1,1) 2 (4,2) 1 (3,4) 5 (4,5) -1 (5,5) -5 将X转化为稀疏存储方式。 X=[2,0,0,0,0; 0,0,0,0,0; 0,0,0,5,0; 0,1,0,0,-1; 0,0,0,0,-5]; A=sparse(X) 其他调用格式： 30

? sparse(m,n)：生成一个m×n的所有元素都是0的稀疏矩阵。 ? sparse(u,v,S)：u,v,S是3个等长的向量

? S是要建立的稀疏矩阵的非0元素。 ? u(i)、v(i)分别是S(i)的行和列下标。

? 该函数建立一个max(u)行、max(v)列并以S为稀疏元素的稀疏矩阵。

和稀疏矩阵操作有关的函数：

? [u,v,S]=find(A)：返回矩阵A中非0元素的下标和元素。这里产生的u,v,S可作为sparse(u,v,S)的参数。

? full(A)：返回和稀疏存储矩阵A对应的完全存储方式矩阵。

2. 产生稀疏存储矩阵

B=spconvert(A)

其中A为一个m×3或m×4的矩阵，其每行表示一个非0元素，m是非0元素的个数，A每个元素的意义是：

(i,1) 第i个非0元素所在的行。 (i,2) 第i个非0元素所在的列。 (i,3) 第i个非0元素值的实部。

(i,4) 第i个非0元素值的虚部，若为实矩阵，则无第四列。

例6.42（产生稀疏存储矩阵）p185

根据表示稀疏矩阵的矩阵A，产生一个稀疏存储方式矩阵B（按列）。

?2?3?A??4??5??6A=[2,2,1; 3,1,-1; 4,3,3; 5,3,8; 6,6,12]; B=spconvert(A) 213361??1??3? ?8?12??B = (3,1) -1 (2,2) 1 (4,3) 3 (5,3) 8 (6,6) 12 3. 带状稀疏存储矩阵

A=spdiags(B,d,m,n)

? m,n为原带状矩阵的行数与列数。 ? B为r×p阶矩阵。 ? r=min(m,n)。

? p为原带状矩阵所有非零对角线的条数。

? B的第i列为原带状矩阵的第i条非零对角线。 ? d为具有p个元素的列向量，它的第i个元素为该带状矩阵的第i条对角线的位置k。 非零对角线的取值方法：

31

? 若对角线上元素个数等于r，则取全部元素；

? 若非零对角线上元素个数小于r，则应该用零补足到r个元素。 补零的原则：

? m

? 若位于主对角线的下方第s条对角线，取k=-s； ? 若位于主对角线的上方第s条对角线，取k=s。

k=-2 -1 0 3×4阶矩阵（m

k=-2 1 0 3×3阶矩阵（m≥n）

-1 -2 0 0 a11 0 0 0 0 1 a11 2 3 0 0 0

a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a12 a13 a21 a22 a23 a31 a32 a33 0 0 0

4×3阶矩阵与3×3阶矩阵类似。 其他调用格式：

? [B,d]=spdiags(A)：从原带状矩阵A中提取全部非零对角线元素赋给矩阵B及其这些非零对角线的位置向量d。 ? B=spdiags(A,d)：从原带状矩阵A中提取由向量d所指定的那些非零对角线元素构成矩阵B。

? E=spdiags(B,d,A)：在原带状矩阵A中将由向量d所指定的那些非零对角线元素用矩阵B替代，构成一个新的带状矩阵E。若赋值号左边改为A，则A为经过替换后的新稀疏矩阵。

例（带状稀疏存储矩阵）p187

矩阵的对角线元素及其位置。 % m

-2 -1 0 1 2 3 0 0 11 12 13 14 0 21 22 23 24 0 31 32 33 34 0 0 % m=n A=[11:13;21:23;31:33] [B,d]=spdiags(A); [d';B] A = 11 12 13 21 22 23 31 32 33 ans = -2 -1 0 1 2 31 21 11 0 0 0 32 22 12 0 0 0 33 23 13 A = 11 12 13 21 22 23 31 32 33 41 42 43 ans = -3 -2 -1 0 1 2 41 31 21 11 0 0 0 42 32 22 12 0 0 0 43 33 23 13 % m>n A=[11:13;21:23; 31:33;41:43] [B,d]=spdiags(A); [d';B] 4. 单位矩阵的稀疏存储

单位矩阵只有对角线元素为1，其他元素都为0，是一种具有稀疏特征的矩阵。 ? eye：产生一个完全存储方式的单位矩阵。

? speye(m,n)：返回一个m×n的稀疏存储单位矩阵。

6.7.3 稀疏矩阵应用举例 p187

稀疏存储矩阵只是矩阵的存储方式不同，它的运算规则与普通矩阵是一样的。所以，在运算过程中，稀疏存储矩阵可直接参与运算。

当参与运算的对象不全是稀疏存储矩阵时，所得结果一般是完全存储形式。

例（稀疏存储矩阵的运算）p187

>> A=[0,0,3; 0,5,0; 0,0,9] A = 0 0 3 0 5 0 0 0 9 >> B=sparse(A) B = (2,2) 5 33

(1,3) 3 (3,3) 9 >> B*B ans = (2,2) 25 (1,3) 27 (3,3) 81 >> A*A ans = 0 0 27 0 25 0 0 0 81 >> rand(3)*B ans = 0 4.7858 4.1716 0 2.4269 4.2687 0 4.0014 11.1534 例6.43（带状稀疏存储矩阵的运算）p188

求三对角线性代数方程组的解。

?23??x1??0??141??x??3????2????164??x3???2? ??????262???x4??1??11?????5???x5???>> B=[ 1,2,0; 1,4,3; 2,6,1; 1,6,4; 0,1,2]; >> d=[-1;0;1]; >> A=spdiags(B,d,5,5) A = (1,1) 2 (2,1) 1 (1,2) 3 (2,2) 4 (3,2) 1 (2,3) 1 (3,3) 6 (4,3) 2 (3,4) 4 (4,4) 6 34

(5,4) 1 (4,5) 2 (5,5) 1 >> b=[0;3;2;1;5]; >> x=(inv(A)*b)' x = -0.1667 0.1111 2.7222 -3.6111 8.6111 >> A=full(A) A = 2 3 0 0 0 1 4 1 0 0 0 1 6 4 0 0 0 2 6 2 0 0 0 1 1 >> x=(inv(A)*b)' x = -0.1667 0.1111 2.7222 -3.6111 8.6111 表 稀疏矩阵函数p183~189

函 数 sparse find full spconvert spdiags eye speye

说 明 完全存储方式转化为稀疏存储方式 返回矩阵中非0元素的下标和元素 返回和稀疏存储矩阵对应的完全存储方式矩阵 产生稀疏存储矩阵 产生带状稀疏存储矩阵 产生完全存储方式的单位矩阵 返回稀疏存储的单位矩阵 35

本文来源：https://www.bwwdw.com/article/35qo.html