现代科学运算-matlab语言与应用
东北大学 http://www.icourses.cn/home/ (薛定宇)
《高等应用数学问题的MATLAB求解(第四版)》
代码可在matlab r2016b 运行。
03 微积分问题的计算机求解
03.01 极限计算
单变量函数的极限
极限的定义:
L=limx→x0f(x)
L
=
lim
x
→
x
0
f
(
x
)
matlab函数: L = limit(fun, x,
x0
x
0
)
左右极限:
L1=limx→x−0f(x)L2=limx→x+0f(x)
L
1
=
lim
x
→
x
0
−
f
(
x
)
L
2
=
lim
x
→
x
0
+
f
(
x
)
matlab函数: L = limit(fun, x,
x0
x
0
, ‘left’ or ‘right’)
例3-1 求极限
limx→0sinxx
lim
x
→
0
sin
x
x
三步:声明符号变量、描述函数、直接求解
% matlab 求极限
syms x; f = sin(x)/x; limit(f, x, 0)
例3-2 求极限 limx→∞x(1+ax)xsinbx lim x → ∞ x ( 1 + a x ) x sin b x
syms x a b; % 变量之间用空格,不能用,
f = x*(1+a/x)^x*sin(b/x);
limit(f, x, inf)
例3-3 求单边极限 limx→0+ex3−11−cosx−sinx−−−−−−−√ lim x → 0 + e x 3 − 1 1 − cos x − sin x
syms x;
limit((exp(x^3)-1)/(1-cos(sqrt(x-sin(x)))), x, 0, 'right')
极限及拓展知识
x = -0.1: 0.001: 0.1;
y = (exp(x .^3) - 1) ./ (1-cos(sqrt(x - sin(x))));
plot(x, y, '-', [0], [12], 'o')
拓展:复数的余弦 – Euler公式
cosjα=(eα+e−α)/2
c
o
s
j
α
=
(
e
α
+
e
−
α
)
/
2
例3-4 正切函数的单边极限
函数 tan t
求关于
π/2
π
/
2
点处的左右极限
syms t;
f = tan(t);
L1 = limit(f, t, pi/2, 'left'),
L2 = limit(f, t, pi/2, 'right')
% 间断点
例3-5 序列的极限
给定序列:
limn→inftyn2−−√3sinn!n+1
lim
n
→
i
n
f
t
y
n
2
3
sin
n
!
n
+
1
syms n;
f = n^(2/3)*sin(factorial(n))/(n+1);
F = limit(f, n, inf)
例3-6 求序列函数极限
limn→∞narctan(1n(x2+1)+x)tann(π4+x2n)
lim
n
→
∞
n
arctan
(
1
n
(
x
2
+
1
)
+
x
)
tan
n
(
π
4
+
x
2
n
)
syms x, n;
f = n*atan(1/(n*(x^2+1)+x))*tan(pi/4+x/2/n)^n;
limit(f, n, inf)
03.02 区间极限与多变量极限
例3-7 求极限 limx→∞xn,limn→∞xn lim x → ∞ x n , lim n → ∞ x n
syms x n real;
f = x^n;
L1 = limit(f, x, inf), L2 = limit(f, n, inf)
MuPAD的底层limit函数
底层limit函数可以求解区间极限问题
L = feval(symengine, ‘limit’, f, ‘x=infinity’, ‘Intervals’)
例3-8 a,b > 0,求 limx→∞asin(8x2)+bcos(2x−2) lim x → ∞ a sin ( 8 x 2 ) + b cos ( 2 x − 2 )
syms a b positive, syms x;
f = a*sin(8*x^2)+b*cos(2*x-2);
L=feval(symengine, 'limit', f, 'x=infinity', 'Intervals')
分段函数的解析描述
% 分段函数的解析描述
function f=piecewise(varargin), str=[];
try
for i=1:2:length(varargin),
str = [str, '[', varargin{i}, ',', varargin[i+1], '],'];
end
catch
error('Input arguments should be given in pairs.'),
end
f = feval(symengine, 'piecewise', str(1:end-1));
% 调用格式
% f = piecewise(var1, var2, ...)
例3-9 分段函数绘图
分段函数
y={1.1sign(x),|x|>1.1x,|x|≤1.1
y
=
{
1.1
s
i
g
n
(
x
)
,
|
x
|
>
1.1
x
,
|
x
|
≤
1.1
f = piecewise('abs(x)>1.1', ...
'1.1*sign(x)', 'abs(x) <= 1.1', 'x');
syms x;
x0 = -3: 0.01: 3;
f1 = subs(f, x, x0);
plot(x0, f1)
多变量函数的极限
函数f(x, y)的累极限
L1=limx→x0[limy→y0f(x,y)],L2=limy→y0[limx→x0f(x,y)]
L
1
=
lim
x
→
x
0
[
lim
y
→
y
0
f
(
x
,
y
)
]
,
L
2
=
lim
y
→
y
0
[
lim
x
→
x
0
f
(
x
,
y
)
]
累极限求解
L1
L
1
= limit(limit(f, y,
y0
y
0
), x,
x0
x
0
)
L2
L
2
= limit(limit(f, x,
x0
x
0
), y,
y0
y
0
)
积分顺序不同,结果可能不同
例3-10 二元函数的累极限
limy→∞[limx→1/y√e−1/(y2+x2)sin2xx2(1+1y2)x+a2y2]
lim
y
→
∞
[
lim
x
→
1
/
y
e
−
1
/
(
y
2
+
x
2
)
sin
2
x
x
2
(
1
+
1
y
2
)
x
+
a
2
y
2
]
syms x a;
syms y positive;
f = exp(-1/(y^2+x^2))*sin(x)^2/x^2 ...
*(1+1/y^2)^(x+a^2*y^2);
L = limit(limit(f, x, 1/sqrt(y)), y, inf)
多元函数重极限
数学表达式
L=limx→x0y→y0f(x,y)
L
=
lim
x
→
x
0
y
→
y
0
f
(
x
,
y
)
不易求解
理论上,沿所有方向均得出相同极限才可
不能用累极限方法求解
累极限存在单不相等,没有重极限
有时累极限存在且相等,但无重极限
例3-11 计算重极限
limx→∞y→∞(xyx2+y2)x2
lim
x
→
∞
y
→
∞
(
x
y
x
2
+
y
2
)
x
2
试图按4个不同方向和速度求
syms x y;
f = (x*y/(x^2+y^2))^(x^2);
L1 = limit(limit(f, x, inf), y, inf),
L2 = limit(limit(f, y, inf), x, inf)
L3 = limit(limit(f, x, y^2), y, inf),
L4 = limit(limit(f, y, x^2), x, inf)
例3-12 判断重极限是否存在
limx→0y→0xyx2+y2
lim
x
→
0
y
→
0
x
y
x
2
+
y
2
沿y = rx 方向趋近
syms r x y;
f = x*y/(x^2+y^2);
L = limit(subs(f, y, r*x), x, 0)
% L = r/(r^2 + 1),r不同,极限不同,重极限不存在
03.03 单变量函数求导
函数的导数和高阶导数
函数求导的三个步骤
df(x)dx=limh→0f(t+h)−f(t)h
d
f
(
x
)
d
x
=
lim
h
→
0
f
(
t
+
h
)
−
f
(
t
)
h
声明变量、输入函数表达式、调用diff()函数直接求导
函数语法
df(x)dx
d
f
(
x
)
d
x
: y = diff(fun, x)
或
dnf(x)dxn
d
n
f
(
x
)
d
x
n
: y = diff(fun, x, n)
自变量为唯一符号变量,可以省去x
例3-13 一阶导数及其曲线
给定函数
f(x)=sinxx2+4x+3
f
(
x
)
=
sin
x
x
2
+
4
x
+
3
, 求
d4f(x)dx4
d
4
f
(
x
)
d
x
4
% 求1阶导数和4阶导数
syms x;
f = sin(x)/(x^2 + 4*x + 3);
f1 = diff(f)
f4 = diff(f, x, 4)
ezplot(f, [0, 5]),
hold on; ezplot(f1, [0, 5]),
hold on; ezplot(f4, [0, 5])
hold off;
% 化简
ff1 = collect(simplify(f4), sin(x))
ff4 = collect(simplify(f4), cos(x))
% 函数diff()的高效率
tic, diff(f, x, 50); toc
例3-14 复合泛函求导
已知函数
F(t)=t2sintf(t)
F
(
t
)
=
t
2
sin
t
f
(
t
)
推导其3阶导函数公式
当
f(t)=e−t
f
(
t
)
=
e
−
t
时,求3阶导数
难点:如何定义f(t) – syms f(t)
syms t f(t);
G = simplify(diff(t^2*sin(t)*f, t, 3))
% f(t) = e^{-t}
simplify(subs(G, f, exp(-t))),
% 验证
simplify(diff(t^2*sin(t)*exp(-t), 3) -ans)
例3-15 矩阵函数求导
矩阵函数
H(x)=[4sin5x3x2+4x+1e−4x24x2+2−−−−−−√]
H
(
x
)
=
[
4
sin
5
x
e
−
4
x
2
3
x
2
+
4
x
+
1
4
x
2
+
2
]
直接求解,对每个矩阵元素直接求导
syms x;
H=[4*sin(5*x), exp(-4*x^2);
3*x^2+4*x+1, sqrt(4*x^2+2)];
H13 = diff(H, x, 3)
参数方程的导数
参数方程:
y=f(t),x=g(t)
y
=
f
(
t
)
,
x
=
g
(
t
)
由递推公式求:
dnydxn
d
n
y
d
x
n
dydx=f′(t)g′(t)
d
y
d
x
=
f
′
(
t
)
g
′
(
t
)
d2ydx2=ddt(f′(t)g′(t))1g′(t)=ddt(dydx)1g′(t)
d
2
y
d
x
2
=
d
d
t
(
f
′
(
t
)
g
′
(
t
)
)
1
g
′
(
t
)
=
d
d
t
(
d
y
d
x
)
1
g
′
(
t
)
⋮
⋮
dnydxn=ddt(dn−1ydxn−1)1g′(t)
d
n
y
d
x
n
=
d
d
t
(
d
n
−
1
y
d
x
n
−
1
)
1
g
′
(
t
)
参数方程求导的递归实现
% 参数方程求导的递归实现函数
function result = paradiff(y, x, t, n)
if mod(n, 1) ~= 0
error('n should positive integer, please correct'),
else
if n == 1
result = diff(y, t)/diff(x, t);
else
result = diff(paradiff(y, x, t, n-1), t)/diff(x, t);
end;
end;
end
% y1 = paradiff(y, x, t, n)
例3-16 参数方程求导
已知参数方程
y=sint(t+1)3,x=cost(t+1)3
y
=
sin
t
(
t
+
1
)
3
,
x
=
cos
t
(
t
+
1
)
3
, 求
d3ydx3
d
3
y
d
x
3
。
syms t;
y = sin(t)/(t+1)^3;
x = cos(t)/(t+1)^3;
f = paradiff(y, x, t, 3);
[n, d] = numden(f); % 提取分子、分母
F = simplify(n)/simplify(d) %分子、分母分别化简,再相除
03.04 偏导数计算
多元函数偏导数
双变量函数
f(x,y)
f
(
x
,
y
)
求高阶偏导数:
∂m+n∂xm∂ynf(x,y)
∂
m
+
n
∂
x
m
∂
y
n
f
(
x
,
y
)
matlab语法:
f = diff(diff(fun, x, m), y, n)
或 f = diff(diff(fun, y, n), x, m)
例3-17 一阶偏导数
二元函数偏导数图形表示
z=f(x,y)=(x2−2x)e−x2−y2−xy
z
=
f
(
x
,
y
)
=
(
x
2
−
2
x
)
e
−
x
2
−
y
2
−
x
y
求两个偏导数:
∂z∂x,∂z∂y
∂
z
∂
x
,
∂
z
∂
y
syms x y;
z = (x^2 - 2*x)*exp(-x^2-y^2-x*y);
zx = simplify(diff(z, x, 1)),
zy = diff(z, y, 1)
梯度函数图形表示
绘制三维曲面
[x0, y0] = meshgrid(-3: 0.2: 2, -2: .2: 2);
z0 = double(subs(z, {x, y}, {x0, y0}));
surf(x0, y0, z0), zlim([-0.7 1.5])
引力线(负梯度)
contour(x0, y0, z0, 30), hold on;
zx0 = subs(zx, {x, y}, {x0, y0});
zy0 = subs(zy, {x, y}, {x0, y0});
quiver(x0, y0, -zx0, -zy0)
例3-18 三元函数求偏导
三元函数
f(x,y,z)=sin(x2y)e−x2y−z2
f
(
x
,
y
,
z
)
=
sin
(
x
2
y
)
e
−
x
2
y
−
z
2
求偏导数:
∂4f(x,y,z)∂x2∂y∂z
∂
4
f
(
x
,
y
,
z
)
∂
x
2
∂
y
∂
z
syms x y z;
f = sin(x^2*y)*exp(-x^2*y-z^2);
df = diff(diff(diff(f, x, 2), y), z);
df = simplify(df)
隐函数偏导数
隐函数:
f(x1,x2,⋯,xn)=0
f
(
x
1
,
x
2
,
⋯
,
x
n
)
=
0
一阶偏导数:
∂xi∂xj=−∂∂xjf(x1,x2,⋯,xn)∂∂xif(x1,x2,⋯,xn)
∂
x
i
∂
x
j
=
−
∂
∂
x
j
f
(
x
1
,
x
2
,
⋯
,
x
n
)
∂
∂
x
i
f
(
x
1
,
x
2
,
⋯
,
x
n
)
matlab代码: F = - diff(f,
xj
x
j
) / diff(f,
xi
x
i
)
二元隐函数的高阶偏导数
二元函数 f(x, y)
一阶偏导数:
∂y/∂x=F1(x,y)
∂
y
/
∂
x
=
F
1
(
x
,
y
)
二阶偏导数:
F2(x,y)=∂2y∂x2=∂F1(x,y)∂x+∂F1(x,y)∂yF1(x,y)
F
2
(
x
,
y
)
=
∂
2
y
∂
x
2
=
∂
F
1
(
x
,
y
)
∂
x
+
∂
F
1
(
x
,
y
)
∂
y
F
1
(
x
,
y
)
高阶偏导数:
Fn(x,y)=∂ny∂xn=∂Fn−1(x,y)∂x+∂Fn−1(x,y)∂yF1(x,y)
F
n
(
x
,
y
)
=
∂
n
y
∂
x
n
=
∂
F
n
−
1
(
x
,
y
)
∂
x
+
∂
F
n
−
1
(
x
,
y
)
∂
y
F
1
(
x
,
y
)
高阶偏导数求解
递推公式:
Fn(x,y)=∂ny∂xn=∂Fn−1(x,y)∂x+∂Fn−1(x,y)∂yF1(x,y)
F
n
(
x
,
y
)
=
∂
n
y
∂
x
n
=
∂
F
n
−
1
(
x
,
y
)
∂
x
+
∂
F
n
−
1
(
x
,
y
)
∂
y
F
1
(
x
,
y
)
matlab求解:
f1
f
1
= impldiff(f, x, y, n)
% 求高阶偏导数
function dy = impldiff(f, x, y, n)
if mod(n, 1) ~= 0
error('n should positive integer, please correct')
else
F1 = -simplify(diff(f, x)/diff(f, y));
dy = F1;
for i=2:n
dy = simplify(diff(dy, x) + diff(dy, y)*F1);
end;
end;
end
例3-19 隐函数求导
二元隐函数
z=f(x,y)=(x2−2x)e−x2−y2−xy
z
=
f
(
x
,
y
)
=
(
x
2
−
2
x
)
e
−
x
2
−
y
2
−
x
y
% 一阶偏导数
syms x y;
f = (x^2-2*x)*exp(-x^2-y^2-x*y);
F1 = impldiff(f, x, y, 1)
% 二阶偏导数
F2 = impldiff(f, x, y, 2)
% 三阶偏导数
F3 = impldiff(f, x, y, 3)
[n, d] = numden(F3), simplify(n)
例3-20 隐函数求导与简化
隐函数:
x2+xy+y2=3
x
2
+
x
y
+
y
2
=
3
各阶导数:
f = x^2 + x*y + y^2 - 3;
F1 = impldiff(f, x, y, 1)
f2 = impldiff(f, x, y, 2);
F2 = subs(f2, x^2+x*y+y^2, 3)
f3 = impldiff(f, x, y, 3);
F3 = subs(f3, x^2+x*y+y^2, 3)
f4 = impldiff(f, x, y, 4);
F4 = subs(f4, x^2+x*y+y^2, 3)
多元函数的Jacobi矩阵
多元函数
⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪y1=f1(x1,x2,⋯,xn)y2=f2(x1,x2,⋯,xn)⋮⋮ym=fm(x1,x2,⋯,xn)
{
y
1
=
f
1
(
x
1
,
x
2
,
⋯
,
x
n
)
y
2
=
f
2
(
x
1
,
x
2
,
⋯
,
x
n
)
⋮
⋮
y
m
=
f
m
(
x
1
,
x
2
,
⋯
,
x
n
)
Jacobi矩阵
J=⎡⎣⎢⎢⎢⎢⎢∂y1/∂x1∂y2/∂x1⋮∂ym/∂x1∂y1/∂x2∂y2/∂x2⋮∂ym/∂x2⋯⋯⋱⋯∂y1/∂xn∂y2/∂xn⋮∂ym/∂xm⎤⎦⎥⎥⎥⎥⎥
J
=
[
∂
y
1
/
∂
x
1
∂
y
1
/
∂
x
2
⋯
∂
y
1
/
∂
x
n
∂
y
2
/
∂
x
1
∂
y
2
/
∂
x
2
⋯
∂
y
2
/
∂
x
n
⋮
⋮
⋱
⋮
∂
y
m
/
∂
x
1
∂
y
m
/
∂
x
2
⋯
∂
y
m
/
∂
x
m
]
matlab求解: J = jacobian(y, x)
例3-21 Jacobi矩阵
直角坐标和极坐标变化公式
x=rsinθcosϕ,y=rsinθsinϕ,z=rcosθ
x
=
r
sin
θ
cos
ϕ
,
y
=
r
sin
θ
sin
ϕ
,
z
=
r
cos
θ
求 Jacobi矩阵
syms r theta phi;
x = r*sin(theta)*cos(phi);
y = r*sin(theta)*sin(phi);
z = r*cos(theta);
J = jacobian([x;y;z], [r theta phi])
Hessian偏导数矩阵
给定n元函数
f(x1,x2,⋯,xn)
f
(
x
1
,
x
2
,
⋯
,
x
n
)
Hessian矩阵
H=⎡⎣⎢⎢⎢⎢⎢∂2f/∂x21∂2f/∂x2∂x1⋮∂2f/∂xn∂x1∂2f/∂x1∂x2∂2f/∂x22⋮∂2f/∂xn∂x2⋯⋯⋱⋯∂2f/∂x1∂xn∂2f/∂x2∂xn⋮∂2f/∂x2n⎤⎦⎥⎥⎥⎥⎥
H
=
[
∂
2
f
/
∂
x
1
2
∂
2
f
/
∂
x
1
∂
x
2
⋯
∂
2
f
/
∂
x
1
∂
x
n
∂
2
f
/
∂
x
2
∂
x
1
∂
2
f
/
∂
x
2
2
⋯
∂
2
f
/
∂
x
2
∂
x
n
⋮
⋮
⋱
⋮
∂
2
f
/
∂
x
n
∂
x
1
∂
2
f
/
∂
x
n
∂
x
2
⋯
∂
2
f
/
∂
x
n
2
]
matlab:H = hessian(f, x)
早期matlab版本: H = jacobian(jacobian(f, x), x)
例3-22 Hessian矩阵
试求二元函数的Hessian矩阵
z=f(x,y)=(x2−2x)e−x2−y2−xy
z
=
f
(
x
,
y
)
=
(
x
2
−
2
x
)
e
−
x
2
−
y
2
−
x
y
syms x y;
f = (x^2-2*x)*exp(-x^2-y^2-x*y);
H = simplify(hessian(f, [x, y]))
% 另一种方法
X = [x, y];
H1 = simplify(jacobian(jacobian(f, X), X))
03.05 积分运算
* 求不定积分*
不定积分:
∫f(x)dx,∫⋯∫f(x)dxn
∫
f
(
x
)
d
x
,
∫
⋯
∫
f
(
x
)
d
x
n
matlab函数:
积分: F = int(fun, x)
多重积分,嵌套调用;更多重应该用循环
F = int(… int(fun, x))
积分得出的是原函数,结果再加C
例3-23 先求导再积分还原
函数:
f(x)=sinxx2+4x+3
f
(
x
)
=
sin
x
x
2
+
4
x
+
3
求其一阶导数,再积分
syms x;
y = sin(x) / (x^2 + 4*x + 3);
y1 = diff(y);
y0 = int(y1)
% 求其4阶导数,再积分,检验结果
y4 = diff(y, 4);
y0 = int(int(int(int(y4))));
simplify(y0)
例3-26 积分公式证明
试证明:
∫x3cos2axdx=x48+(x34a−3x8a3)sin2ax+(3x28a2−316a4)cos2ax+C
∫
x
3
cos
2
a
x
d
x
=
x
4
8
+
(
x
3
4
a
−
3
x
8
a
3
)
sin
2
a
x
+
(
3
x
2
8
a
2
−
3
16
a
4
)
cos
2
a
x
+
C
% 对等号左侧进行化简
syms a x;
f = simplify(int(x^3*cos(a*x)^2, x));
% 输入右侧,比较并化简
f1 = x^4/8 + (x^3/(4*a) - 3*x/(8*a^3))*sin(2*a*x) + (3*x^2/(8*a^2) - 3/(16*a^4))*cos(2*a*x);
simplify(f-f1)
例3-27 不可积函数的不定积分计算
考虑如下两个不可积问题
f(x)=e−x2/2,g(x)=xsin(ax4)ex2/2
f
(
x
)
=
e
−
x
2
/
2
,
g
(
x
)
=
x
sin
(
a
x
4
)
e
x
2
/
2
matlab求解:
syms x;
I = int(exp(-x^2/2))
结果:
I=erf(2–√x)2π−−√
I
=
e
r
f
(
2
x
)
2
π
特殊函数:
erf(x)=2π−−√∫x0e−t2dt
e
r
f
(
x
)
=
2
π
∫
0
x
e
−
t
2
d
t
含参数的不可积函数
求解
∫g(x)dx
∫
g
(
x
)
d
x
, 其中
g(x)=xsin(ax4)ex2/2
g
(
x
)
=
x
sin
(
a
x
4
)
e
x
2
/
2
syms a x;
int(x*sin(a*x^4)*exp(x^2/2))
% 并不是所有的积分都能被计算出,因为原始函数不一定存在
% ans = int(x*sin(a*x^4)*exp(x^2/2), x)
定积分与无穷积分计算
定积分与无穷积分
∫baf(x)dx∫∞af(x)dx
∫
a
b
f
(
x
)
d
x
∫
a
∞
f
(
x
)
d
x
语句格式:
I = int(fun, x, a, b)
I = int(fun, x, a, inf)
有时需要Newton-Leibniz公式 - 原函数 F(x)
∫baf(x)dx=F(b)−F(a)
∫
a
b
f
(
x
)
d
x
=
F
(
b
)
−
F
(
a
)
例3-28 定积分求解
函数
f(x)=e−x2/2
f
(
x
)
=
e
−
x
2
/
2
求a = 0, b = 1.5 或
∞
∞
时的定积分值
syms x;
I1 = int(exp(-x^2/2), x, 0, 1.5)
vpa(I1, 70),
I2 = int(exp(-x^2/2), x, 0, inf)
例3-30 反常积分计算
数学问题
∫2e11x1−ln2x−−−−−−−√dx
∫
1
2
e
1
x
1
−
ln
2
x
d
x
syms x;
f = 1/(x*sqrt(1 - log(x)^2));
I = int(f, x, 1, 2*exp(sym(1)))
多重积分问题
多重积分
I=∫xMxm∫yM(x)ym(x)∫zM(x,y)zm(x,y)f(x,y,z)dzdydx
I
=
∫
x
m
x
M
∫
y
m
(
x
)
y
M
(
x
)
∫
z
m
(
x
,
y
)
z
M
(
x
,
y
)
f
(
x
,
y
,
z
)
d
z
d
y
d
x
利用int()函数直接求解
注意:某些问题需要根据实际情况先选择积分顺序,可积的部分作为内积分,然后再处理外积分,否则会的不出解析解
例3-32 三元函数积分
已知三元函数
F(x,y,z)=−4ze−x2y−z2[cosx2y−10yx2cosx2y+4x4y2sinx2y+4x4y2cosx2y−sinx2y]
F
(
x
,
y
,
z
)
=
−
4
z
e
−
x
2
y
−
z
2
[
cos
x
2
y
−
10
y
x
2
cos
x
2
y
+
4
x
4
y
2
sin
x
2
y
+
4
x
4
y
2
cos
x
2
y
−
sin
x
2
y
]
试求
∫⋯∫F(x,y,z)dx2dydz
∫
⋯
∫
F
(
x
,
y
,
z
)
d
x
2
d
y
d
z
syms x y z;
f0 = -4*z*exp(-x^2*y-z^2)*(cos(x^2*y) - 10*cos(x^2*y)*y*x^2 ...
+4*sin(x^2*y)*x^4*y^2 + 4*cos(x^2*y)*x^4*y^2 - sin(x^2*y));
f1 = int(f0, z);
f1 = int(f1, y);
f1 = int(f1, x);
f1 = simplify(int(f1, x))
% f1 = sin(x^2*y)*exp(- y*x^2 - z^2)
例3-33 多重积分
求解积分问题
∫20∫π0∫π04xze−x2y−z2dzdydx
∫
0
2
∫
0
π
∫
0
π
4
x
z
e
−
x
2
y
−
z
2
d
z
d
y
d
x
syms x y z;
I = int(int(int(4*x*z*exp(-x^2*y-z^2), z, 0, pi), y, 0, pi), x, 0, 2)
% I=-(2*exp(-pi^2) - 2)*(eulergamma/2 + log(2) + log(pi)/2 - ei(-4*pi)/2)
注意:
eulergamma 为euler常数
γ
γ
Ei(n,z)=∫∞1e−ztt−ndt
E
i
(
n
,
z
)
=
∫
1
∞
e
−
z
t
t
−
n
d
t
03-06 Fourier级数逼近
Fourier级数展开
给定周期函数
f(x),xin[−L,L],T=2L
f
(
x
)
,
x
i
n
[
−
L
,
L
]
,
T
=
2
L
满足:
f(x)=f(x+kT)
f
(
x
)
=
f
(
x
+
k
T
)
一般
(a,b)
(
a
,
b
)
区间,变量替换
L=b−a2,x1=x+L+a
L
=
b
−
a
2
,
x
1
=
x
+
L
+
a
Fourier级数展开
f(x)=a02+∑n=1∞(ancosnπLx+bnsinnπLx)
f
(
x
)
=
a
0
2
+
∑
n
=
1
∞
(
a
n
cos
n
π
L
x
+
b
n
sin
n
π
L
x
)
系数:
⎧⎩⎨⎪⎪⎪⎪an=1L∫L−Lf(x)cosnπxLdx,n=0,1,2,⋯bn=1L∫L−Lf(x)sinnπxLdx,n=1,2,3,⋯
{
a
n
=
1
L
∫
−
L
L
f
(
x
)
cos
n
π
x
L
d
x
,
n
=
0
,
1
,
2
,
⋯
b
n
=
1
L
∫
−
L
L
f
(
x
)
sin
n
π
x
L
d
x
,
n
=
1
,
2
,
3
,
⋯
Fourier级数展开函数编写
求解 [F, A, B] = fseries(f, x, p, a, b)
% 求Fourier级数展开的matlab代码
function [F, A, B] = fseries(f, x, varargin)
[p, a, b] = default_vals({6, -pi, pi}, varargin{:});
L = (b-a)/2;
if a+b
f = subs(f, x, x+L+a);
end;
A = int(f, x, -L, L)/L;
B = [];
F = A/2;
for n = 1:p
an = int(f*cos(n*pi*x/L), x, -L, L)/L;
A = [A, an];
bn = int(f*sin(n*pi*x/L), x, -L, L)/L;
B = [B, bn];
F = F + an*cos(n*pi*x/L) + bn*sin(n*pi*x/L);
end;
if a+b
F = subs(F, x, x-L-a);
end;
end
function varargout = default_vals(vals, varargin)
if nargout ~= length(vals)
error('number of arguments mismatch');
else
nn = length(varargin) + 1;
varargout = varargin;
for i = nn:nargout
varargout{i} = vals{i};
end;
end;
end
例3-33 函数的Fourier近似
函数
y=x(x−π)(x−2π),x∈(0,2π)
y
=
x
(
x
−
π
)
(
x
−
2
π
)
,
x
∈
(
0
,
2
π
)
求它的Fourier级数展开
syms x;
f = x*(x-pi)*(x - 2*pi);
[F, A, B] = fseries(f, x, 12, 0, 2*pi);
F
%比较结果
ezplot(f, [0, 2*pi]),
hold on,
ezplot(F, [0, 2*pi]), hold off;
展开结果和效果
更大的区间
x∈(−π,3π)
x
∈
(
−
π
,
3
π
)
ezplot(f, [-pi, 3*pi]),
hold on,
ezplot(F, [-pi, 3*pi]),
hold off;
例3-34 分段函数的近似
给定函数
y={1x≥0−1x<0x∈(−π,π)
y
=
{
1
x
≥
0
−
1
x
<
0
x
∈
(
−
π
,
π
)
原函数可以表示成
f(x)=|x|x
f
(
x
)
=
|
x
|
x
syms x;
f = abs(x)/x;
xx = [-pi: pi/200: pi];
xx = xx(xx ~= 0);
xx = sort([xx, -eps, eps]);
yy = subs(f, x, xx);
plot(xx, yy),
hold on;
for n = 1: 20
[f1, a, b] = fseries(f, x, n);
y1 = subs(f1, x, xx);
plot(xx, y1)
end;
展开效果与函数
前14项的Fourier级数展开
[f1, a, b] = fseries(f, x, 14); f1
数学形式:
f(x)≈4sinxπ+4sin3x3π+4sin5x5π+4sin7x7π+4sin9x9π+4sin11x11π+4sin13x13π
f
(
x
)
≈
4
s
i
n
x
π
+
4
sin
3
x
3
π
+
4
sin
5
x
5
π
+
4
sin
7
x
7
π
+
4
sin
9
x
9
π
+
4
sin
11
x
11
π
+
4
sin
13
x
13
π
一般形式:
f(x)=4π∑k=1∞sin(2k−1)x2k−1
f
(
x
)
=
4
π
∑
k
=
1
∞
sin
(
2
k
−
1
)
x
2
k
−
1
周期函数假设
在区间
[−2π,2π]
[
−
2
π
,
2
π
]
进行拟合效果比较
xx = [-2*pi: pi/200:2*pi];
xx = xx(xx ~= 0);
xx = sort([xx, -eps, eps]);
yy = subs(f, x, xx);
plot(xx, yy),
hold on;
y1 = subs(f1, x, xx);
plot(xx, y1),
hold off;
开始假设: f(t)=f(kT+t) f ( t ) = f ( k T + t )
03.07 Taylor级数逼近
单变量函数的Taylor幂级数展开
数学表示
在 x = 0 点附近的Taylor 幂级数
f(x)=a1+a2x+a3x2+⋯+akxk−1+o(xk)
f
(
x
)
=
a
1
+
a
2
x
+
a
3
x
2
+
⋯
+
a
k
x
k
−
1
+
o
(
x
k
)
其中
ai=1i!limx→0di−1dxi−1f(x),i=1,2,3,⋯
a
i
=
1
i
!
lim
x
→
0
d
i
−
1
d
x
i
−
1
f
(
x
)
,
i
=
1
,
2
,
3
,
⋯
matlab格式: F = taylor(fun, x, ‘Order’, k)
早期版本: F = taylor(fun, x, k)
关于 x = a 的 Taylor 展开
关于 x = a 点的 Taylor 展开
f(x)=b1+b2(x−a)+b3(x−a)2+⋯+bk(x−a)k−1+o([(x−a)k]
f
(
x
)
=
b
1
+
b
2
(
x
−
a
)
+
b
3
(
x
−
a
)
2
+
⋯
+
b
k
(
x
−
a
)
k
−
1
+
o
(
[
(
x
−
a
)
k
]
其中
bi=1i!limx→adi=1dxi−1f(x),i=1,2,3,⋯
b
i
=
1
i
!
lim
x
→
a
d
i
=
1
d
x
i
−
1
f
(
x
)
,
i
=
1
,
2
,
3
,
⋯
matlab格式: F = taylor(fun, x, a, ‘Order’, k)
早期版本: F = taylor(fun, x, a, k)
例3-35 函数的Taylor展开
函数
f(x)=sinxx2+4x+3
f
(
x
)
=
sin
x
x
2
+
4
x
+
3
在 x = 0, x = 2 和 x = a 处的 Taylor 展开.
% 在 x = 0 的 Taylor 展开
syms x;
f = sin(x)/(x^2 + 4*x + 3);
y = taylor(f, x, 'Order', 9)
% 在x = 2 的 Taylor 展开
y = taylor(f, x, 2, 'Order', 9)
% 在x = a 的 Taylor 展开
syms a;
y = taylor(f, x, a, 'Order', 9)
近似效果
% 在 x = 0 的 Taylor 展开
syms x;
f = sin(x)/(x^2 + 4*x + 3);
y = taylor(f, x, 'Order', 9)
% 在区间[-1, 1]内绘图
ezplot(f, [-1, 1]),
hold on;
h = ezplot(y, [-1, 1]);
set(h, 'Color', 'r'),
hold off;
% 更小的区间
figure;
ezplot(f, [-0.6, 0.6]),
hold on;
h = ezplot(y, [-0.6, 0.6]);
set(h, 'Color', 'r');
hold off;
例3-36 正弦函数的幂级数逼近
原函数
y(t)=sin(t)
y
(
t
)
=
sin
(
t
)
syms x;
y = sin(x);
ezplot(y),
hold on;
for n = [6: 2: 16]
p = taylor(y, x, 'Order', n),
ezplot(p);
end;
hold off;
多变量函数的Taylor级数展开
多元函数
f(x1,x2,⋯,xn)
f
(
x
1
,
x
2
,
⋯
,
x
n
)
Taylor幂级数展开
f(x)=f(a)+[(x1−a1)∂∂x1+⋯+(xn−an)∂∂xn]f(x)|x=a+
f
(
x
)
=
f
(
a
)
+
[
(
x
1
−
a
1
)
∂
∂
x
1
+
⋯
+
(
x
n
−
a
n
)
∂
∂
x
n
]
f
(
x
)
|
x
=
a
+
12![(x1−a1)∂∂x1+⋯+(xn−an)∂∂xn]2f(x)|x=a+⋯+
1
2
!
[
(
x
1
−
a
1
)
∂
∂
x
1
+
⋯
+
(
x
n
−
a
n
)
∂
∂
x
n
]
2
f
(
x
)
|
x
=
a
+
⋯
+
1k![(x1−a1)∂∂x1+⋯+(xn−an)∂∂xn]kf(x)|x=a+⋯
1
k
!
[
(
x
1
−
a
1
)
∂
∂
x
1
+
⋯
+
(
x
n
−
a
n
)
∂
∂
x
n
]
k
f
(
x
)
|
x
=
a
+
⋯
F = taylor(f, [
x1,x2,⋯,xn
x
1
,
x
2
,
⋯
,
x
n
], [
a1,x2,⋯,an
a
1
,
x
2
,
⋯
,
a
n
], ‘Order’, k)
例3-37 二元函数Taylor展开
函数
z=f(x,y)=(x2−2x)e−x2−y2−xy
z
=
f
(
x
,
y
)
=
(
x
2
−
2
x
)
e
−
x
2
−
y
2
−
x
y
% 在原点[0, 0]展开Taylor级数
syms x y;
f = (x^2-2*x)*exp(-x^2-y^2-x*y);
F = taylor(f, [x, y], [0, 0], 'Order', 9);
F1 = collect(F, x)
% 关于在(1, a)点展开
syms a;
F = taylor(f, [x,y], [1, a], 'Order', 9);
F2 = simplify(F)
03.08 级数求和与序列求积
级数求和
已知通项的有穷或无穷级数的和
数学表示:
S=∑k=k0knfk
S
=
∑
k
=
k
0
k
n
f
k
matlab语句: S = symsum(
fk,k,k0,kn
f
k
,
k
,
k
0
,
k
n
)
例子3-38 由古老的故事引出级数求和问题
计算
S=20+2122+⋯+262+263=∑i=0632i
S
=
2
0
+
2
1
2
2
+
⋯
+
2
6
2
+
2
6
3
=
∑
i
=
0
63
2
i
format long;
sum(2 .^ [0:63])
% 使用symsum()
syms k;
symsum(2^k, 0, 63)
% 更多项的扩展
sum(sym(2) .^ [0:200])
例3-39 级数求和
求无穷级数的和
S=11×4+14×7+⋯+1(3n−2)(3n+1)+⋯
S
=
1
1
×
4
+
1
4
×
7
+
⋯
+
1
(
3
n
−
2
)
(
3
n
+
1
)
+
⋯
syms n;
s = symsum(1/((3*n-2)*(3*n+1)), n, 1, inf)
% 使用数值方法
m = 1:10000000;
s1 = sum(1 ./((3*m-2) .* (3*m+1)));
format long;
s1
例3-40 级数函数的求和
试求解含有变量x的无穷级数
J=2∑n=0∞1(2n+1)(2x+1)2n+1
J
=
2
∑
n
=
0
∞
1
(
2
n
+
1
)
(
2
x
+
1
)
2
n
+
1
符号运算方法
早期版本不能给出级数的收敛区域
syms x n;
s1 = symsum(2/((2*n+1)*(2*x+1)^(2*n+1)), n, 0, inf);
simplify(s1)
% piecewise(0 < x | x < -1, 2*atanh(1/(2*x + 1)))
例3-41 综合问题
试求级数与极限综合问题
limn→∞[(1+12+13+14+⋯+1n)−lnn]
lim
n
→
∞
[
(
1
+
1
2
+
1
3
+
1
4
+
⋯
+
1
n
)
−
ln
n
]
syms m n;
res = limit(symsum(1/m, m, 1, n) - log(n), n, inf)
vpa(res)
% 注意:该问题不能先求解无穷级数的和,在减去 ln n,都是无穷大
% ans = eulergamma
例3-42 级数的极限
S=limn→∞[(1+1n2sinπn2+⋯+(1+n−1n2)sin(n−1)πn2]
S
=
lim
n
→
∞
[
(
1
+
1
n
2
sin
π
n
2
+
⋯
+
(
1
+
n
−
1
n
2
)
sin
(
n
−
1
)
π
n
2
]
通项提取:
(1+kn2)sinkπn2,k=1,2,⋯,n−1
(
1
+
k
n
2
)
sin
k
π
n
2
,
k
=
1
,
2
,
⋯
,
n
−
1
syms n k;
S = simplify(limit(symsum((1+k/n^2)*sin(k*pi/n^2), k, 1, n-1), n, inf))
序列求积问题
序列求积运算:
P=∏n=abf(n)
P
=
∏
n
=
a
b
f
(
n
)
matlab语句:
新符号运算工具箱: P = symprod(
fn
f
n
, n, a, b)
早期版本: P = maple(‘product’, fun, ‘n=a..b’)
例3-43 序列求积
计算序列乘积:
Pn=∏k=1n(1+1k3)
P
n
=
∏
k
=
1
n
(
1
+
1
k
3
)
% 序列求积
syms k n;
P1 = symprod(1 + 1/k^3, k, 1, n);
P1 = simplify(P1)
P2 = symprod(1 + 1/k^3, 1, inf);
P2 = simplify(P2)
例3-44 序列求积
S=1−12+1×32×4−1×3×52×4×6+1×3×5×72×4×6×8−1×3×5×7×92×4×6×8×10+⋯
S
=
1
−
1
2
+
1
×
3
2
×
4
−
1
×
3
×
5
2
×
4
×
6
+
1
×
3
×
5
×
7
2
×
4
×
6
×
8
−
1
×
3
×
5
×
7
×
9
2
×
4
×
6
×
8
×
10
+
⋯
通项:
(−1)n∏k=1n2k−12k,k=1,2,⋯,∞
(
−
1
)
n
∏
k
=
1
n
2
k
−
1
2
k
,
k
=
1
,
2
,
⋯
,
∞
syms k n;
S = 1 + symsum((-1)^n*symprod((2*k-1)/(2*k), k, 1, n), n, 1, inf)
% 2^(1/2)/2
03.09 曲线积分与曲面积分
第一类曲线积分
第一类曲线积分:
I1=∫lf(x,y,z)ds
I
1
=
∫
l
f
(
x
,
y
,
z
)
d
s
曲线l满足参数方程 x = x(t), y = y(t), z = z(t)
ds=(dxdt)2+(dydt)2+(dzdt)2−−−−−−−−−−−−−−−−−−−√dt
d
s
=
(
d
x
d
t
)
2
+
(
d
y
d
t
)
2
+
(
d
z
d
t
)
2
d
t
转换成定积分问题
I=∫tMtmf[x(t),y(t),z(t)]x2t+y2t+z2t−−−−−−−−−−√dt
I
=
∫
t
m
t
M
f
[
x
(
t
)
,
y
(
t
)
,
z
(
t
)
]
x
t
2
+
y
t
2
+
z
t
2
d
t
曲线积分函数
% 通用曲线积分函数
function I = path_integral(F, vars, t, a, b)
if length(F) == 1
I = int(F*sqrt(sum(diff(vars, t) .^2)), t, a, b);
else
F = F(:) .';
vars = vars(:);
I = int(F*diff(vars, t), t, a, b);
end;
end
调用格式:
I = path_integral(f, [x, y], t,
tm,tM
t
m
,
t
M
)
I = path_integral(f, [x, y, z], t,
tm,tM
t
m
,
t
M
)
例3-49 曲线积分求解
计算
∫lz2x2+y2ds,l
∫
l
z
2
x
2
+
y
2
d
s
,
l
是如下定义的螺线
⎧⎩⎨⎪⎪x=acost,y=asint,z=at,0≤t≤2π,a>0
{
x
=
a
cos
t
,
y
=
a
sin
t
,
z
=
a
t
,
0
≤
t
≤
2
π
,
a
>
0
公式:
I=∫tMtmf[x(t),y(t),z(t)]x2t+y2t+z2t−−−−−−−−−−√dt
I
=
∫
t
m
t
M
f
[
x
(
t
)
,
y
(
t
)
,
z
(
t
)
]
x
t
2
+
y
t
2
+
z
t
2
d
t
% 计算曲线积分
syms t;
syms a positive;
x = a*cos(t);
y = a* sin(t);
z = a*t;
f = z^2/(x^2+y^2);
I = path_integral(f, [x, y, z], t, 0, 2*pi)
% (8*2^(1/2)*a*pi^3)/3
第二类曲线积分
第二类曲线积分
I2=∫lf⃗ (x,y,z)⋅ds⃗
I
2
=
∫
l
f
→
(
x
,
y
,
z
)
⋅
d
s
→
其中
f⃗ (x,y,z)=[P(x,y,z),Q(x,y,z),R(x,y,z)]
f
→
(
x
,
y
,
z
)
=
[
P
(
x
,
y
,
z
)
,
Q
(
x
,
y
,
z
)
,
R
(
x
,
y
,
z
)
]
并且:
ds⃗ =[dxdt,dydt,dzdt]Tdt
d
s
→
=
[
d
x
d
t
,
d
y
d
t
,
d
z
d
t
]
T
d
t
上式化为:
∫ba[P(x,y,z),Q(x,y,z),R(x,y,z)][dxdt,dydt,dzdt]Tdt
∫
a
b
[
P
(
x
,
y
,
z
)
,
Q
(
x
,
y
,
z
)
,
R
(
x
,
y
,
z
)
]
[
d
x
d
t
,
d
y
d
t
,
d
z
d
t
]
T
d
t
matlab语句:
I = path_integral([P, Q], [x, y], t, a, b)
I = path_integral([P, Q, R], [x, y, z], t, a, b)
I = path_integral(F, v, t, a, b)
例3-51 第二类曲线积分
曲线积分:
∫lx+yx2+y2dx−x−yx2+y2dy
∫
l
x
+
y
x
2
+
y
2
d
x
−
x
−
y
x
2
+
y
2
d
y
l为正向圆周x2+y2=a2
l
为
正
向
圆
周
x
2
+
y
2
=
a
2
正向圆周的参数函数描述
x=acost,y=asint,(0≤t≤2π)
x
=
a
cos
t
,
y
=
a
sin
t
,
(
0
≤
t
≤
2
π
)
syms t;
syms a positive;
x = a * cos(t);
y = a * sin(t);
F = [(x+y)/(x^2+y^2), -(x-y)/(x^2+y^2)];
I = path_integral(F, [x, y], t, 2*pi, 0)
曲面积分
% 曲面积分函数
function I=surf_integral(f,xx,uu,um,vm)
if length(f)==1 % scalar integrand
if length(xx)==1 % surface described by an explicit function
I=int(int(f*sqrt(1+diff(xx,uu(1))^2+diff(xx,uu(2))^2),...
uu(2),um(1),um(2)),uu(1),vm(1),vm(2));
else % surface described by an implicit function
xx=[xx(:).' 1]; x=xx(1); y=xx(2); z=xx(3); u=uu(1); v=uu(2);
E=diff(x,u)^2+diff(y,u)^2+diff(z,u)^2;
F=diff(x,u)*diff(x,v)+diff(y,u)*diff(y,v)+diff(z,u)*diff(z,v);
G=diff(x,v)^2+diff(y,v)^2+diff(z,v)^2;
I=int(int(f*sqrt(E*G-F^2),u,um(1),um(2)),v,vm(1),vm(2));
end
else % vector integrand
if length(xx)==1 % surface described by an explicit function
syms x y z; ua=sqrt(1+diff(xx,x)^2+diff(xx,y)^2);
cA=-diff(xx,x)/ua; cB=-diff(xx,y)/ua; cC=1/ua;
I=surf_integral(f(:).'*[cA; cB; cC],xx,uu,um,vm);
else, x=xx(1); y=xx(2); z=xx(3); u=uu(1); v=uu(2);
A=diff(y,u)*diff(z,v)-diff(z,u)*diff(y,v);
B=diff(z,u)*diff(x,v)-diff(x,u)*diff(z,v);
C=diff(x,u)*diff(y,v)-diff(y,u)*diff(x,v); F=A*f(1)+B*f(2)+C*f(3);
I=int(int(F,uu(1),um(1),um(2)),uu(2),vm(1),vm(2));
end;
end
第一类曲面积分
I = surf_integral(
f,z,[x,y],[ym,yM],[xm,xM]
f
,
z
,
[
x
,
y
]
,
[
y
m
,
y
M
]
,
[
x
m
,
x
M
]
)
I = surf_integral(
f,[x,y,z],[u,v],[um,uM],[vm,vM]
f
,
[
x
,
y
,
z
]
,
[
u
,
v
]
,
[
u
m
,
u
M
]
,
[
v
m
,
v
M
]
)
第二类曲面积分
I = surf_integral(
[P,Q,R],z,[u,v],[um,uM],[vm,vM]
[
P
,
Q
,
R
]
,
z
,
[
u
,
v
]
,
[
u
m
,
u
M
]
,
[
v
m
,
v
M
]
)
I = surf_integral(
[P,Q,R],[x,y,z],[u,v],[um,uM],[vm,vM]
[
P
,
Q
,
R
]
,
[
x
,
y
,
z
]
,
[
u
,
v
]
,
[
u
m
,
u
M
]
,
[
v
m
,
v
M
]
)
例3-54 曲面积分
曲面积分
∬(x2y+zy2)ds
∬
(
x
2
y
+
z
y
2
)
d
s
积分曲面
x=ucosv,y=usinv,z=v,0≤u≤a,0≤v≤2π
x
=
u
cos
v
,
y
=
u
sin
v
,
z
=
v
,
0
≤
u
≤
a
,
0
≤
v
≤
2
π
syms u v;
syms a positive;
x = u * cos(v);
y = u * sin(v);
z = v;
f = x^2*y + z*y^2;
I = surf_integral(f, [x, y, z], [u, v], [0, a], [0, 2*pi])
03.10 数值微分
一阶数值微分算法
导数定义:
y′(t)=limΔt→0y(t+Δt)−y(t)Δt
y
′
(
t
)
=
lim
Δ
t
→
0
y
(
t
+
Δ
t
)
−
y
(
t
)
Δ
t
前向差分公式:
y′i≈ΔyiΔt=yi+1−yiΔt
y
i
′
≈
Δ
y
i
Δ
t
=
y
i
+
1
−
y
i
Δ
t
后向差分公式:
y′i≈ΔyiΔt=yi−yi−1Δt
y
i
′
≈
Δ
y
i
Δ
t
=
y
i
−
y
i
−
1
Δ
t
算法精度:
o(Δt)
o
(
Δ
t
)
中心差分算法高精度公式
y′i≈−yi+2+8yi+1−8yi−1+yi−212Δt
y
i
′
≈
−
y
i
+
2
+
8
y
i
+
1
−
8
y
i
−
1
+
y
i
−
2
12
Δ
t
y′′i≈−yi+2+16yi+1−30yi+16yi−1−yi−212Δt2
y
i
′
′
≈
−
y
i
+
2
+
16
y
i
+
1
−
30
y
i
+
16
y
i
−
1
−
y
i
−
2
12
Δ
t
2
y′′′i≈−yi+3+8yi+2−13yi+1+13yi−1−8yi−2+yi−38Δt3
y
i
′
′
′
≈
−
y
i
+
3
+
8
y
i
+
2
−
13
y
i
+
1
+
13
y
i
−
1
−
8
y
i
−
2
+
y
i
−
3
8
Δ
t
3
y(4)i≈−yi+3+12yi+2−39yi+1+56yi−39yi−1+12yi−2−yi−36Δt4
y
i
(
4
)
≈
−
y
i
+
3
+
12
y
i
+
2
−
39
y
i
+
1
+
56
y
i
−
39
y
i
−
1
+
12
y
i
−
2
−
y
i
−
3
6
Δ
t
4
算法精度
o(Δt4)
o
(
Δ
t
4
)
中心差分算法
% 中心差分算法
function [dy,dx]=diff_ctr(y,Dt,n)
y1=[y 0 0 0 0 0 0];
y2=[0 y 0 0 0 0 0];
y3=[0 0 y 0 0 0 0];
y4=[0 0 0 y 0 0 0];
y5=[0 0 0 0 y 0 0];
y6=[0 0 0 0 0 y 0];
y7=[0 0 0 0 0 0 y]; %{\rm\,尽量采用向量化运算,准备若干平移的向量}
switch n %{\rm\,按给定的阶次由开关结构求数值微分}
case 1, dy=(-y1+8*y2-8*y4+y5)/12/Dt;
case 2, dy=(-y1+16*y2-30*y3+16*y4-y5)/12/Dt^2;
case 3, dy=(-y1+8*y2-13*y3+13*y5-8*y6+y7)/8/Dt^3;
case 4, dy=(-y1+12*y2-39*y3+56*y4-39*y5+12*y6-y7)/6/Dt^4;
end
dy=dy(5+2*(n>2):end-4-2*(n>2)); dx=([2:length(dy)+1]+(n>2))*Dt;
end
函数调用: [dy,dx] [ d y , d x ] = diff_ctr( y,Δt,n y , Δ t , n )
例3-56 数值微分逼近
对函数
f(x)=sinxx2+4x+3
f
(
x
)
=
sin
x
x
2
+
4
x
+
3
h = 0.05;
x = 0: h: pi;
syms x1;
y = sin(x1)/(x1^2+4*x1+3);
yy1 = diff(y);
f1 = subs(yy1, x1, x);
yy2 = diff(yy1);
f2 = subs(yy2, x1, x);
yy3 = diff(yy2);
f3 = subs(yy3, x1, x);
yy4 = diff(yy3);
f4 = subs(yy4, x1, x);
与解析解的精确比较
% 比较不同阶导数
y = sin(x1)/(x1^2+4*x1+3);
y = subs(y, x1, x);
[y1, dx1] = diff_ctr(y, h, 1);
subplot(221), plot(x, f1, dx1, y1, ':');
[y2, dx2] = diff_ctr(y, h, 2);
subplot(222), plot(x, f2, dx2, y2, ':');
[y3, dx3] = diff_ctr(y, h, 3);
subplot(223), plot(x, f3, dx3, y3, ':');
[y4, dx4] = diff_ctr(y, h, 4);
subplot(224), plot(x, f4, dx4, y4, ':');
% 分析误差
norm(double((y4-f4(4:60))./f4(4:60)))
二元函数梯度计算
多变量函数
f(x1,x2,⋯,xn)
f
(
x
1
,
x
2
,
⋯
,
x
n
)
梯度是其对各个自变量偏导数构成的向量
v=[∂f∂x1,∂f∂x2,⋯,∂f∂xn]
v
=
[
∂
f
∂
x
1
,
∂
f
∂
x
2
,
⋯
,
∂
f
∂
x
n
]
函数gradient()的调用格式
[fx,fy]=gradient(z)
[
f
x
,
f
y
]
=
g
r
a
d
i
e
n
t
(
z
)
计算梯度
fx=fx/Δx,fy=fy/Δy
f
x
=
f
x
/
Δ
x
,
f
y
=
f
y
/
Δ
y
其中Δx和Δy分别是x和y生成网格的步距
其
中
Δ
x
和
Δ
y
分
别
是
x
和
y
生
成
网
格
的
步
距
例3-57 梯度数值计算
已知
z=f(x,y)=(x2−2x)e−x2−y2−xy
z
=
f
(
x
,
y
)
=
(
x
2
−
2
x
)
e
−
x
2
−
y
2
−
x
y
生成数据
用数值方法求梯度并分析误差
syms x y;
f = (x^2-2*x)*exp(-x^2-y^2-x*y);
[x1, y1] = meshgrid(-3: .2: 3, -2: .2: 2);
z1 = double(subs(f, {x, y}, {x1, y1}));
[zx1, zy1] = gradient(z1);
zx1 = zx1/0.2;
zy1 = zy1/0.2;
contour(x1, y1, z1, 30);
hold on;
quiver(x1, y1, -zx1, -zy1),
hold off;
03.11 单变量函数的数值积分
由给定数据进行梯形求积
梯形近似方法的基本思想
∫baf(x)dx=∑i=1N∫xi+1xif(x)dx=∑i=1NΔfi
∫
a
b
f
(
x
)
d
x
=
∑
i
=
1
N
∫
x
i
x
i
+
1
f
(
x
)
d
x
=
∑
i
=
1
N
Δ
f
i
matlab调用格式: sum((2*y(1:end-1, :) + diff(y)) .* diff(x)) /2
或 S = trapz(x, y)
例3-59 定步长求积分
用定步长方法求解积分
∫3π20cos15xdx
∫
0
3
π
2
cos
15
x
d
x
并比较不同步长下的结果
首先绘图
% 在求解区域内被积函数有很强的震荡
x = [0:0.01:3*pi/2, 3*pi/2];
y = cos(15*x);
plot(x, y)
不同步距积分结果
syms x, A = int(cos(15*x), 0, 3*pi/2);
h0 = [0.1, 0.01, 0.001, 0.0001, 0.00001, 0.000001];
v = [];
for h=h0
x = [0:h:3*pi/2, 3*pi/2];
y = cos(15*x);
I = trapz(x, y);
v = [v; h, I, 1/15-I];
end;
v
数值积分求解函数
调用格式
新版本数值积分
I = integral(f, a, b, property pairs)
早期版本
求定积分
[y,k]=quad(fun,a,b)
[
y
,
k
]
=
q
u
a
d
(
f
u
n
,
a
,
b
)
[y,k]=quad(fun,a,b,ε)
[
y
,
k
]
=
q
u
a
d
(
f
u
n
,
a
,
b
,
ε
)
其它数值定积分函数:quadl(), quadgk(), quadv()
例3-60 被积函数描述
计算积分: erf(x) =
2π−−√∫1.50e−t2dt
2
π
∫
0
1.5
e
−
t
2
d
t
方法1:一般函数方法(点运算)
function y = c3ffun(x)
y = 2/sqrt(pi)*exp(-x.^2);
方法2:inline函数方法(不建议使用)
f = inline('2/sqrt(pi)*exp(-x.^2)', 'x');
方法3:匿名函数(matlab 7.0+)
f = @(x)2/sqrt(pi)*exp(-x.^2);
数值求解
y = integral(f, 0, 1.5)
解析解:运用符号工具箱
syms x;
y0 = vpa(int(2/sqrt(pi)*exp(-x^2), 0, 1.5), 60)
默认选项下数值解函数 integral()即可保证高精度的数值解
f = @(x)2/sqrt(pi)*exp(-x.^2);
y = integral(f, 0, 1.5)
syms x;
y0 = vpa(int(2/sqrt(pi)*exp(-x^2), 0, 1.5), 60)
例3-61 分段函数积分
给定如下分段函数
f(x)=⎧⎩⎨⎪⎪ex2,0≤x≤2804−sin(16πx),2<x≤4
f
(
x
)
=
{
e
x
2
,
0
≤
x
≤
2
80
4
−
sin
(
16
π
x
)
,
2
<
x
≤
4
计算积分值:
I=∫40f(x)dx
I
=
∫
0
4
f
(
x
)
d
x
定积分示意图:
x = [0:0.01:2, 2+eps: 0.01: 4, 4];
y = exp(x.^2).*(x <= 2) + ...
80./(4-sin(16*pi*x)).*(x>2);
y(end) = 0;
x = [eps, x];
y = [0,y]
fill(x, y, 'g')
求解与验证
% 求解语句
f = @(x)exp(x.^2).*(x<=2) + 80 ./ ...
(4-sin(16*pi*x)).*(x>2);
I = integral(f, 0, 4)
% 提高精度(检验)
I2 = integral(f, 0, 4, 'RelTol', 1e-20)
%解析解
syms x;
f = piecewise(x<=2, exp(x^2), ...
x>2, 80/(4-sin(16*pi*x)));
I3 = vpa(int(f, x, 0, 4))
例3-62 与梯形法比较
重新计算积分
∫3π/20cos(15x)dx
∫
0
3
π
/
2
cos
(
15
x
)
d
x
f = @(x) cos(15*x);
tic,
S = integral(f, 0, 3*pi/2, 'RelTol', 1e-20),
toc
结论:和梯形法相比, 速度今读都明显提高
理论值: 1/15
例3-64 大范围积分
积分问题:
∫1000cos(15x)dx
∫
0
100
c
o
s
(
15
x
)
d
x
早期版本的几个积分函数不适用
% 数值解
f = @(x)cos(15*x);
I1 = integral(f, 0, 100, 'RelTol', 1e-20)
% 解析解
syms x;
I = int(cos(15*x), x, 0, 100),
vpa(I)
广义数值积分问题求解
Integral()函数可以直接用于广义积分运算
函数调用格式与前面介绍的完全一致
I = integral(fun, a, b)
早期版本, 采用Gauss-Kronrod算法
[y, k] = quadgk(fun, a, b,
ε
ε
)
例3-66 广义积分的数值计算
数值计算:
∫∞0e−x2dx
∫
0
∞
e
−
x
2
d
x
数值解与解析解
解析解不可求,可以引出特殊函数erf()
数值解可以直接得出数值
f = @(x)exp(-x.^2);
I = integral(f, 0, inf, 'RelTol', 1e-20)
syms x;
I1 = int(exp(-x^2), 0, inf),
vpa(I1)
例3-67 含参数函数积分
积分函数: 绘制I(a)曲线
I(a)=∫∞0eax2sin(a2x)dx
I
(
a
)
=
∫
0
∞
e
a
x
2
sin
(
a
2
x
)
d
x
新版本支持向量参数a
早期版本需要对a循环
% 含参数函数积分
a = 0: 0.1: 4;
f = @(x) exp(-a*x.^2).*sin(a.^2*x);
I = integral(f, 0, inf, 'RelTol', 1e-20, 'ArrayValued', true);
plot(a, I)
03.12 双重数值积分
双重积分问题的数值解
数学标准型
I=∫xMxm∫yM(x)ym(x)f(x,y)dydx
I
=
∫
x
m
x
M
∫
y
m
(
x
)
y
M
(
x
)
f
(
x
,
y
)
d
y
d
x
注意积分顺序,先y后x
I = integral2(f,
xm,xM,ym,yM
x
m
,
x
M
,
y
m
,
y
M
, par pairs)
ym与yM还可以为函数句柄
y
m
与
y
M
还
可
以
为
函
数
句
柄
例3-68 双重积分
试求出双重定积分
J=∫1−1∫2−2e−x2/2sin(x2+y)dxdy
J
=
∫
−
1
1
∫
−
2
2
e
−
x
2
/
2
sin
(
x
2
+
y
)
d
x
d
y
积分区域为矩形区域, 对应好边界即可
f = @(x, y) exp(-x.^2/2).*sin(x.^2+y);
J = integral2(f, -2, 2, -1, 1, 'RelTol', 1e-20)
积分曲面绘制
仿照积分曲线绘制的函数,
[x, y, F] = intfunc2(f,
xm,xM,ym,yM
x
m
,
x
M
,
y
m
,
y
M
, n, m)
% 积分曲面
function [yv, xv, F] = intfunc2(f, xm, xM, varargin)
[ym, yM, n, m] = default_vals({xm,xM,50,50},varargin{:});
xv = linspace(xm, xM, n);
yv = linspace(ym, yM, m);
d = yv(2) - yv(1);
[x y] = meshgrid(xv, yv);
F = zeros(n,m);
for i = 2:n
for j = 2:m
F(i,j) = integral2(f, xv(1), xv(i), yv(1), yv(j), 'RelTol', 1e-20);
end;
end
例3-69 二元函数积分曲面
绘制积分函数曲面
J=∫1−1∫2−2e−x2/2sin(x2+y)dxdy
J
=
∫
−
1
1
∫
−
2
2
e
−
x
2
/
2
sin
(
x
2
+
y
)
d
x
d
y
积分区域为矩形区域, 对应好边界即可
f = @(x, y) exp(-x.^2/2).*sin(x.^2+y);
[x, y, z] = intfunc2(f, -2, 2, -1, 1);
surf(x, y, z);
I = z(end, end)
例3-70 双重积分计算
是求出双重定积分
J=∫1−12∫1−x2/2√−1−x2/2√e−x22sin(x2+y)dydx
J
=
∫
−
1
2
1
∫
−
1
−
x
2
/
2
1
−
x
2
/
2
e
−
x
2
2
sin
(
x
2
+
y
)
d
y
d
x
先y后x,与标准型一致, 可以直接求解
fh = @(x) sqrt(1-x.^2/2);
fl = @(x) -sqrt(1-x.^2/2);
f = @(x, y) exp(-x.^2/2).*sin(x.^2+y);
y = integral2(f, -1/2, 1, fl, fh)
先x后y的数值积分计算
与标准型积分顺序相反的二元积分
I=∫yMym∫xM(y)xm(y)f(x,y)dxdy
I
=
∫
y
m
y
M
∫
x
m
(
y
)
x
M
(
y
)
f
(
x
,
y
)
d
x
d
y
被积函数描述时替换次序, 令
x^=y,y^=x
x
^
=
y
,
y
^
=
x
则
I=∫x^Mx^m∫y^M(x^)y^m(x^)f(y^,x^)dy^dx^
I
=
∫
x
^
m
x
^
M
∫
y
^
m
(
x
^
)
y
^
M
(
x
^
)
f
(
y
^
,
x
^
)
d
y
^
d
x
^
被积函数: f = @(y, x)
例3-71 解析不可积
先x后y积分
J=∫1−1∫1−y2√−1−y2√e−x2/2sinh(x2+y)dxdy
J
=
∫
−
1
1
∫
−
1
−
y
2
1
−
y
2
e
−
x
2
/
2
sinh
(
x
2
+
y
)
d
x
d
y
解析不可求解,只能求数值解
交换x,y次序
tic,
f = @(y, x)exp(-x.^2/2).*sinh(x.^2+y);
fh = @(y) sqrt(1-y.^2);
fl = @(y)-sqrt(1-y.^2);
I = integral2(f, -1, 1, fl, fh, 'RelTol', 1e-20),
toc;
例3-72 非矩形区域积分区域的积分曲面
数学模型
J=∫1−1∫1−y2√−1−y2√ex2/2sin(x2+y)dxdy
J
=
∫
−
1
1
∫
−
1
−
y
2
1
−
y
2
e
x
2
/
2
s
i
n
(
x
2
+
y
)
d
x
d
y
绘制积分曲面
先计算矩形区域,再删去区域外函数值(NaN)
f = @(x, y)exp(-x.^2/2).*sin(x.^2+y);
fh = @(x)sqrt(1-x.^2/2);
fl = @(x)-sqrt(1-x.^2/2);
[x, y, z] = intfunc2(f, -1/2, 1, -1.2, 1.2);
for i = 1:length(x)
x0 = x(i);
xx = sort([fl(x0), fh(x0)]);
ii = find(y>xx(2) | y < xx(1));
z(ii, i) = NaN;
end;
surf(x, y, z);
03.13 三重与多重数值积分
三重定积分数值求解
长方体区域的三重积分标准型
I=∫xMxm∫yM(x)ym(x)∫zM(x,y)zm(x,y)f(x,y,z)dzdydx
I
=
∫
x
m
x
M
∫
y
m
(
x
)
y
M
(
x
)
∫
z
m
(
x
,
y
)
z
M
(
x
,
y
)
f
(
x
,
y
,
z
)
d
z
d
y
d
x
函数调用:
I=integral3(f,xm,xM,ym,yM,zm,zM,pars)
I
=
i
n
t
e
g
r
a
l
3
(
f
,
x
m
,
x
M
,
y
m
,
y
M
,
z
m
,
z
M
,
p
a
r
s
)
早期版本:
I=triplequad(f,xm,xM,ym,yM,zm,zM,ϵ,@quadl)
I
=
t
r
i
p
l
e
q
u
a
d
(
f
,
x
m
,
x
M
,
y
m
,
y
M
,
z
m
,
z
M
,
ϵ
,
@
q
u
a
d
l
)
例3-73 三重积分计算
用数值方法求三重定积分问题
∫20∫π0∫π04xze−x2y−z2dzdydx
∫
0
2
∫
0
π
∫
0
π
4
x
z
e
−
x
2
y
−
z
2
d
z
d
y
d
x
长方体区域
f = @(x, y, z)4*x.*z.*exp(-x.*x.*y - z.*z);
I = integral3(f, 0, 2, 0, pi, 0, pi, 'RelTol', 1e-20)
多重积分数值求解
NIT工具箱(数值积分工具箱)还可以解决多重超长方体边界定积分问题
I=∫x1Mx1m∫x2Mx2m⋯∫xpMxpmf(x1,x2,⋯,xp)dxp⋯dx2dx1
I
=
∫
x
1
m
x
1
M
∫
x
2
m
x
2
M
⋯
∫
x
p
m
x
p
M
f
(
x
1
,
x
2
,
⋯
,
x
p
)
d
x
p
⋯
d
x
2
d
x
1
调用格式:
I=quadndg(fun,[x1m,x2m,⋯,xpm],[x1M,x2M,⋯,xpM],ϵ)
I
=
q
u
a
d
n
d
g
(
f
u
n
,
[
x
1
m
,
x
2
m
,
⋯
,
x
p
m
]
,
[
x
1
M
,
x
2
M
,
⋯
,
x
p
M
]
,
ϵ
)
例3-74 重新计算三重积分
前面的长方体区域三重积分问题
∫20∫π0∫π04xze−x2y−z2dzdydx
∫
0
2
∫
0
π
∫
0
π
4
x
z
e
−
x
2
y
−
z
2
d
z
d
y
d
x
变量替换:
x1=x,x2=y,x3=z
x
1
=
x
,
x
2
=
y
,
x
3
=
z
被积函数:
f(x)=4x1x3e−x21x2−x23
f
(
x
)
=
4
x
1
x
3
e
−
x
1
2
x
2
−
x
3
2
重新求解(更快):
f = @(x) 4*x(1)*x(3)*exp(-x(1)^2*x(2)-x(3)^2);
tic,
I = quadndg(f, [0 0 0], [2, pi, pi]),
toc
例3-75 五重积分的数值计算
5 重定积分问题
I=∫50∫40∫10∫20∫30v√3w−−√x2y3zdzdydxdwdv
I
=
∫
0
5
∫
0
4
∫
0
1
∫
0
2
∫
0
3
v
3
w
x
2
y
3
z
d
z
d
y
d
x
d
w
d
v
变量替换:
x1=v,x2=w,x3=x,x4=y,x5=z
x
1
=
v
,
x
2
=
w
,
x
3
=
x
,
x
4
=
y
,
x
5
=
z
被积函数:
f(x)=x1−−√3x2−−√x23x34x5
f
(
x
)
=
x
1
3
x
2
x
3
2
x
4
3
x
5
% 五重积分计算
f = @(x)(x(1))^(1/3)*sqrt(x(2))*x(3)^2*x(4)^3*x(5);
I = quadndg(f, [0 0 0 0 0], [5, 4, 1, 2, 3])
五重积分解析解
本例存在解析解
syms x y z w v;
F = v^(1/3)*sqrt(w)*x^2*y^3*z;
I = int(int(int(int(int(F,z,0,3), y, 0, 2), x, 0, 1), w,0, 4), v, 0, 5)
注意:该工具箱单重积分函数duadg()调用格式和quad()一致,其效率也高于duadl(),故在进行数值求积分时建议使用此工具箱。
例3-76 解析不可积的五重积分计算
解析不可积多重积分
I=∫50∫40∫10∫20∫30(e−v√3sinw−−√+e−x2y3z)dzdydxdwdv
I
=
∫
0
5
∫
0
4
∫
0
1
∫
0
2
∫
0
3
(
e
−
v
3
sin
w
+
e
−
x
2
y
3
z
)
d
z
d
y
d
x
d
w
d
v
变量替换:
x1=v,x2=w,xe=x,x4=y,x5=z
x
1
=
v
,
x
2
=
w
,
x
e
=
x
,
x
4
=
y
,
x
5
=
z
被积函数:
f(x)=e−x1√3sinx2−−√+e−x2ex34x5
f
(
x
)
=
e
−
x
1
3
sin
x
2
+
e
−
x
e
2
x
4
3
x
5
f = @(x) exp(-(x(1))^(1/3))*sin(sqrt(x(2)))+exp(-x(3)^2*x(4)^3*x(5));
tic,
I = quadndg(f, [0 0 0 0 0], [5, 4, 1, 2, 3])
toc