第四章数值计算

y0=[1;0];
[t,y]=ode45(@DyDt,tspan,y0);
plot(t,y(:,1)) %x1=x(t)的解
xlabel('t'),title('x(t)') 4
3
2
（4）x(2)为dx/dt的解
dx/dt
1 0
-1
plot(t,y(:,2))
-2
xlabel('t'),ylabel('dx/dt')
©通M信A与T电L子A工B程语学院言 电子教案
源自文库
（2）编写M函数文件 [DyDt.m] function ydot=DyDt(t,x) mu=2; ydot=[x(2);mu*(1-x(1)^2)*x(2)-x(1)];
©通M信A与T电L子A工B程语学院言 电子教案
（3）解算微分方程
tspan=[0,30];
1
fval =
output =
3.1416
iterations: 21
一般解题时后两个输 出参数可以省略！
funcCount: 22 algorithm: [1x46 char]
me©s通Ms信A与aT电gL子Ae工B程: 语学[1院言x电1子1教2 案char]
（2）图形法
xx=-pi/2:pi/200:pi/2;
fun是目标函数名，x1和x2限定自变量的取值范
围。四个输出参数分别表示：极小值点；函数
值；判断是否成功搜到极值点（是否大于0）；
算法。
（2）求多元函数的极值点
[x,fval,exitflag,output]=fminsearch(fun,
x0,option)
求多元函数fun从x0开始的极小值x,该点的函数
©通M信A与T电L子A工B程语学院言 电子教案
4.1.5常微分方程的数值解
Matlab指令格式： [t,y]= ode45(‘F’,tspan,y0) 以上括号中：F是函数文件名；tspan是求数值 解的时间区间；y0是初始条件列向量；t是所求 数值解的自变量列向量,而Y(:,k)就是yk分量的解。 注意：对于高阶微分方程，需要运用变量替换转 换为一阶微分方程，初始条件也作相应的变换。
F= 123 456 789
FX = 111
FX_2 = 222 222 222
FY_2 =
©通M信A与T电L子A工B程语学院言 电子教案
（1）x=eps;
L1 = Ls1 =
L1=(1-cos(2*x))/(x*sin(x)), L2=sin(x)/x
02
L2 = Ls2 =
（2）syms t %可靠的符号计算
f1=(1-cos(2*t))/(t*sin(t));
11
f2=sin(t)/t;
注意：尽量不用
Ls1=limit(f1,t,0)
数值计算求极限
Ls2=limit(f2,t,0)
©通M信A与T电L子A工B程语学院言 电子教案
【例4.1-2】已知x=sin(t)，求该函数在区间
[0, 2] 中的近似导函数。
legend('x(t)','dx/dt')
xlabel('t')
©通M信A与T电L子A工B程语学院言 电子教案
（2）增量取适当
d=pi/100;
t=0:d:2*pi;
x=sin(t);
x_d=sin(t+d);
dxdt_d=(x_d-x)/d;
plot(t,x,'LineWidth',5)
hold on plot(t,dxdt_d) hold off
%小矩形面积之和
s2=d*trapz(y); %折线下面积,可用s_ta=trapz(t,y)代替
s3=d*cumtrapz(y); %折线下累计面积
disp(['sum积分',blanks(3),'trapz积分‘,’ cumtrapz积分’
disp([s1,s2,s3])
©通M信A与T电L子A工B程语学院言 电子教案
…
（3）检查目标函Co数lu值mns 1 through 3 format short e 2.4112e-010 5.7525e+002 2.2967e+003 disp([ff(sx(:,1))C,fof(lusmxn(:4,2)),ff(sx(:,3)),ff(sx(:,4))])
3.3211e+005
增量适当，所求 导函数光滑
legend('x(t)','dx/dt')
xlabel('t') 注意：数值导数的使用应十分谨慎，
自变量增量的取值应适当，符号求导合适。
©通M信A与T电L子A工B程语学院言 电子教案
clf d=pi/100; t=0:d:2*pi; x=sin(t); dxdt_diff=diff(x)/d; dxdt_grad=gradient(x)/d; hold on plot(t,dxdt_grad,'m','LineWidth',8) plot(t(1:end-1),dxdt_diff,'.k','MarkerSize',8)
第四章 数值计算 4.1 数值微积分 4.2 矩阵和代数方程 4.4 多项式运算和卷积
©通M信A与T电L子A工B程语学院言 电子教案
4.1 数值微积分
1.近似数值极限及导数
Matlab数值计算中，没有提供专门的求极限
和求导的指令，但提供了与求导数相关的“求
差分”指令。
（1）DX=diff(x,n,dim)
-3
-4
0
5
10
15
20
25
30
t
©通M信A与T电L子A工B程语学院言 电子教案
（1）高阶变换一阶： 令y1=y, y2=y’,y3=y’’ , 条件为y1(0)=-1, y2(0)=1, y3(0)=0 （2）编写方程函数名为odet3.m（在m文件编辑器里）
function dy=odet3(t,y)
t2=[t,t(end)+d];
y2=[y,nan];
stairs(t2,y2,':k')
hold on
plot(t,y,'r','LineWidth',3)
h=stem(t,y,'LineWidth',2);
set(h(1),'MarkerSize',10)
axis([0,pi/2+d,0,1.5])
yxx=(xx+pi).*exp(abs(sin(xx+pi)));
plot(xx,yxx) xlabel('x'),grid on
13
12
11
10
9
8
7
6
5
4
3
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
x
[xx,yy]=ginput(1)
xx =
%执行后吧十字点对准极值点点击
-6.4885e-007
t=0:pi/100:2*pi; x=sin(t); dt=5*eps; x_eps=sin(t+dt);
增量过小引起有效数 字严重丢失后的毛刺 曲线
dxdt_eps=(x_eps-x)/dt;
plot(t,x,'LineWidth',5)
hold on
plot(t,dxdt_eps)
hold off
或s_n=dblquad(‘x.^y’,0,1,1,2 )
s_n =
0.405466267243508
©通M信A与T电L子A工B程语学院言 电子教案
4.1.4 函数极值的数值求解
（1）求一元函数的最小值点
[x,fval,exitflag,output]=fminbnd(fun,x1,
x2,option)
©通M信A与T电L子A工B程语学院言 电子教案
clear sum求得积分 trapz求得积分 cumtrapz求得积分
d=pi/8; 1.5762 1.3013 t=0:d:pi/2;
0 0.1537 0.4462 0.8450 1.3013
y=0.2+sin(t);
s=sum(y);
s1=d*s;
值为fval.
©通M信A与T电L子A工B程语学院言 电子教案
（1）优化算法
x1=-pi/2; x2=pi/2;
yx=@(x)(x+pi)*exp(abs(sin(x+pi)));
[xn0,fval,exitflag,output]=fminbnd(yx,x1,x2)
xn0 =
exitflag =
-1.2999e-005
786 b= 12 15 15 c=
123
579 ©通M信A与T电L子A工B程语学院言1电2 子教15案 15
3：St=trapz(x,y) 给出采样点（x,y）所连折线下的面积，即 函数y在自变量区间x上的近似积分. 4：Sct=cumtrapz(x,y) 求累计积分 如：s=dt*cumtrapz(Ft); %计算从0开始到每个采样点为止的区间内， Ft曲线下的面积，这个面积是由一个个宽 度为dt的小梯形面积累加而得的。
hold off
shg
注意：以上求积分办
法精度难控制，不常
用。
©通M信A与T电L子A工B程语学院言 电子教案
4.1.3计算精度可控的数值积分
S1=quad(fun,a,b,tol,) Simpson法 s1=quad1(fun,a,b,tol,) lobatto法 功能：a,b为上下限，tol为绝对误差，默认时 积分的绝对精度为10-6
在dim指定的维上，求n阶差分
F= 123
【例】
456
F=[1,2,3;4,5,6;7,8,9]
789
Dx=diff(F)
Dx_2=diff(F,1,2)
%相邻列一阶差分
Dx =
333 ©通M信A与T电L子A工B程语学院言 电子教案
（2）数值梯度 FX =gradient(F,h) 求一元函数的行元素梯度。F为向量时， FX(1)=[F(2)-F(1)]/h FX(end)=[F(end)-F(end-1)]/h, FX(2:end-1)=[(F(3:end)-F(1:end-2))/2]/h h是沿坐标取点的步长，默认步长为1。
©通M信A与T电L子A工B程语学院言 电子教案
4.1.2数值求和与近似数值积分
指令：
1：SX=sum(X)
按列求和，生成一个行向量，各元素是X各
列元素之和。
2：Scs=cumsum(X)
沿列方向求累计和，Scs大小与X的相同，
Sxs(i,k)是X的第k列的 前i个元素之和。
a= 123 456
[例] a=[1 2 3;4 5 6;7 8 6], b=sum(a),c=cumsum(a)
[FX,FY]=gradient(F，h) 求二元函数的梯度，FX,FY 与F的大小相同，分
别表示列、行元素间的“梯度”，h是沿坐标取
点的步长，默认步长为1。
©通M信A与T电L子A工B程语学院言 电子教案
【例】 F=[1,2,3;4,5,6;7,8,9] [FX,FY]=gradient(F) [FX_2,FY_2]=gradient(F,0.5)
Itrapz =
©通M信A与T电L子A工B程语学院言 电子教案
(1)符号计算法
syms x y
s=vpa(int(int(x^y,x,0,1),y,1,2))
s=
.40546510810816438197801311546432
（2）数值积分法
匿名函数解释见P150
format long
s_n=dblquad(@(x,y)x.^y,0,1,1,2)
yy =
©通M信A与T电L子A工B程语学院言 电子教案
sx =
（1）用匿名函数表1.0示00测0 试-0.6函89数7 0.4151 8.0886
ff=@(x)(100*(x(2)1-.x00(100).^-12.)9^1628+(14-.9x6(413))^72.8)0;04 %（注2）意编用写单，纯不形用法xss2fe,求yv.x4a表it1l多1==示2元e-0函10 数极小值点 x0=[-5,-2,2,5;-5,-21,2,5]; %提供4个搜索起点 [sx,sfval,sexit,ssoouutptputu=t]=fminsearch(ff,x0)
©通M信A与T电L子A工B程语学院言 电子教案
方法1：
syms x
Isym=vpa(int(exp(-x^2),x,0,1))
方法2：
format long
d=0.001;x=0:d:1;
Itrapz=d*trapz(exp(-x.*x))
方法3：
Isym =
fx='exp(-x.^2)'; Ic=quad(fx,0,1,1e-8) .7468241328124270253994
S2=dblquad(f,xmin,xmax,ymin,ymax,tol,) 该函数求f(x,y)在[xmin,xmax]×[ymin,ymax]区 域上的二重定积分。参数tol的用法与函数quad 完全相同
S3=triplequad(f,xmin,xmax,ymin,ymax,zmin,
zmax,tol) 该函数求f(x,y)在区域上的三重定积分。参数 tol的用法与函数quad完全相同