现代科学运算-matlab语言与应用
东北大学 http://www.icourses.cn/home/ (薛定宇)
《高等应用数学问题的MATLAB求解(第四版)》
代码可在matlab r2016b 运行。
01 绪论
01.01-02 例子
例1-1 求解: 19931993 1993 1993
%% 1993^{1993}
n = sym(1993)^1993, vpa(n)
例1-2 求函数 f(x)=sinxx2+4x+3的4阶导数d4f(x)dx4 f ( x ) = sin x x 2 + 4 x + 3 的 4 阶 导 数 d 4 f ( x ) d x 4
syms x; f = sin(x)/(x^2 + 4*x + 3);
diff(f, 4)
%%f的50阶导数,显示运行时间
tic, F = diff(f, 50); toc
例1-3 高阶矩阵的行列式如何计算?
80X80的Hilbert矩阵。
H=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢112⋮1n1213⋮1n+11314⋮1n+2⋯⋯⋱⋯1n1n+1⋮12n−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
H
=
[
1
1
2
1
3
⋯
1
n
1
2
1
3
1
4
⋯
1
n
+
1
⋮
⋮
⋮
⋱
⋮
1
n
1
n
+
1
1
n
+
2
⋯
1
2
n
−
1
]
%% 80 * 80 Hilbert matrix
H = hilb(80); H = sym(H); det(H)
例1-4 微分方程求解
⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪y′′+μ(y2−1)y′+y=0,μ=1000,y(0)=−1,y′(0)=1
{
y
′
′
+
μ
(
y
2
−
1
)
y
′
+
y
=
0
,
μ
=
1000
,
y
(
0
)
=
−
1
,
y
′
(
0
)
=
1
mu = 1000;
f = @(t, x)[x(2); -mu*(x(1)^2 - 1) * x(2) - x(1)];
[t, x] = ode15s(f, [0, 3000], [-1, 1]); plot(t, x)
延迟微分方程:
dy(t)dt=−0.1y(t)+0.2y(t−30)1+y10(t−30)
d
y
(
t
)
d
t
=
−
0.1
y
(
t
)
+
0.2
y
(
t
−
30
)
1
+
y
10
(
t
−
30
)
分数阶数微分方程:
3D0.9y(t)3+0.2D0.8y(t)+0.9D0.2y(t)+|2D0.7y(t)|1.5+43y(t)=5sin(10t)
3
D
0.9
y
(
t
)
3
+
0.2
D
0.8
y
(
t
)
+
0.9
D
0.2
y
(
t
)
+
|
2
D
0.7
y
(
t
)
|
1.5
+
4
3
y
(
t
)
=
5
sin
(
10
t
)
例1-5 线性规划问题的求解
min(−2x1−x2−4x3−3x4−x5)
min
(
−
2
x
1
−
x
2
−
4
x
3
−
3
x
4
−
x
5
)
xs.t.⎧⎩⎨⎪⎪2x2+x3+4x4+2x5≤543x1+4x2+5x3−x4−x5≤62x1,x2≥0,x3≥3.32,x4≥0.678,x5≥2.57
x
s
.
t
.
{
2
x
2
+
x
3
+
4
x
4
+
2
x
5
≤
54
3
x
1
+
4
x
2
+
5
x
3
−
x
4
−
x
5
≤
62
x
1
,
x
2
≥
0
,
x
3
≥
3.32
,
x
4
≥
0.678
,
x
5
≥
2.57
f = -[2 1 4 3 1]'; A = [0 2 1 4 2; 3 4 5 -1 -1];
B = [54; 62]; Ae = []; Be = [];
xm = [0, 0, 3.32, 0.678, 2.57];
x = linprog(f, A, B, Ae, Be, xm)
01.03 解析解与数值解
解析解:通过严格推导得到的数学问题的精确解,又称闭式解。
数值解:用数值计算的方法得到的解,通常是问题的近似解。
解析解不存在:
2π−−√∫a0e−x2dx
2
π
∫
0
a
e
−
x
2
d
x
数学家的解决方案 – 发明特殊函数:erf(a)
工程技术人员采用近似解。
01.04 计算机数学语言的发展概述
1978 Cleve Moler, New Mexico University, matlab(Matrix Laboratory)。
三个代表: matlab、Mathematica, Maple。
matlab 仿真: Simulink。
01.05 matlab 解决了计算机语言默认数据类型存储结构有限的问题。
不容易发生因为存储类型引起的溢出错误。
例1-9 Fibonacci数列
a = [1 1];
for i = 3:100, a(i) = a(i-1) + a(i-2); end;
a(end)
%% 用符号类型更精确
a = sym([1 1]);
for i = 3: 300, a(i) = a(i-1) + a(i-2); end;
a(end)
例1-10 矩阵乘法通用子程序
cij=∑k=1paikbkj,i=1,⋯,n,j=1,⋯,m
c
i
j
=
∑
k
=
1
p
a
i
k
b
k
j
,
i
=
1
,
⋯
,
n
,
j
=
1
,
⋯
,
m
%% matlab 会判断是否满足矩阵A的列数等于矩阵B的行数
C = A * B
01.06 三步求解方法
问题是什么、如何描述、求解。
例1-11a 求函数序列的极限
数学问题:
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)
例1-11b 线性规划问题求解
min−2x1−x2−4x3−3x4−x5
min
−
2
x
1
−
x
2
−
4
x
3
−
3
x
4
−
x
5
xs.t.⎧⎩⎨⎪⎪2x2+x3+4x4+2x+5≤543x1+4x2+5x3−x4−x5≤62x1,x2≥0,x3≥3.32,x4≥0.678,x5≥2.57
x
s
.
t
.
{
2
x
2
+
x
3
+
4
x
4
+
2
x
+
5
≤
54
3
x
1
+
4
x
2
+
5
x
3
−
x
4
−
x
5
≤
62
x
1
,
x
2
≥
0
,
x
3
≥
3.32
,
x
4
≥
0.678
,
x
5
≥
2.57
第一步:是什么– 物理含义, x是一个决策向量,求在约束条件下目标函数的最小值。
第二步:如何描述
clear; P.f=[-2 -1 -4 -3 -1];
P.Aineq = [0 2 1 4 2; 3 4 5 -1 -1];
P.Bineq = [54 62]; P.lb = [0; 0; 3.32; 0.678; 2.57];
P.solver='linprog'; P.options = optimoptions('linprog');
第二步:求解
x = linprog(P)
线性混合整数规划
第一步:
物理含义:决策变量的一部分必须是整数,其它任意
如要求1,2,4变量必须为整数。
第二步:描述
clear; P.f=[-2 -1 -4 -3 -1];P.intcon=[1,2,4];
P.Aineq = [0 2 1 4 2; 3 4 5 -1 -1];
P.Bineq = [54 62]; P.lb = [0; 0; 3.32; 0.678; 2.57];
P.solver='intlinprog'; P.options = optimoptions('intlinprog');
第三步:求解
intlinprog(P)
例1-12 神经网络拟合
已知样本点数据向量x, y
x = 0: .5 : 10;
%% 向量的乘、除、幂运算需要在运算符前加点号, [.*] 不能有空格
y = 0.12 * exp(-0.213 * x) + 0.54 * exp(-0.17 * x) .* sin(1.23 * x);
用人工神经网络拟合数据
如何描述:
net = fitnet(5); net.trainParam.epochs = 1000;
求解:
[net, b] = train(net, x, y); plotperform(b)
效果如何:
x0 = [0: 0.1: 10];
y0 = 0.12 * exp(-0.213 * x0) + 0.54 * exp(-0.17 * x0) .* sin(1.23 * x0);
y1 = net(x0); plot(x, y, 'o', x0, y0, x0, y1, ':');