实验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

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