现代科学运算-matlab语言与应用

$\qquad \qquad$ 东北大学 http://www.icourses.cn/home/ (薛定宇)
《高等应用数学问题的MATLAB求解（第四版）》
代码可在matlab r2016b 运行。

01 绪论

01.01-02 例子

例1-1 求解:$1993^{1993}$

%% 1993^{1993}
n = sym(1993)^1993, vpa(n)

例1-2 求函数 $f(x) = \dfrac{\sin x}{x^2 + 4x + 3}的4阶导数 \dfrac{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 = \begin{bmatrix}1 & \dfrac{1}{2} & \dfrac{1}{3} & \cdots & \dfrac{1}{n} \\ \dfrac{1}{2} & \dfrac{1}{3} & \dfrac{1}{4} & \cdots & \dfrac{1}{n+1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \dfrac{1}{n} & \dfrac{1}{n+1} & \dfrac{1}{n+2} & \cdots & \dfrac{1}{2n-1} \end{bmatrix}$

%% 80 * 80 Hilbert matrix
H = hilb(80); H = sym(H); det(H)

例1-4 微分方程求解
$\left \{ \begin{array}{l}y^{\prime \prime} + \mu(y^2 -1) y^{\prime} + y = 0, \\ \mu = 1000, \\ y(0) = -1, \\ y^{\prime}(0) = 1 \end{array} \right.$

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)

延迟微分方程:
$\dfrac{dy(t)}{dt} = -0.1 y(t) + 0.2 \dfrac{y(t - 30)}{1 + y^{10}(t - 30)}$
分数阶数微分方程:
$\dfrac{3 \mathscr{D}^{0.9} y(t)}{3 + 0.2 \mathscr{D}^{0.8} y(t) + 0.9 \mathscr{D}^{0.2} y(t)} + |2 \mathscr{D}^{0.7} y(t)|^{1.5} + \dfrac{4}{3} y(t) = 5 \sin(10t)$

例1-5 线性规划问题的求解
$\min \qquad (-2x_1 - x_2 - 4x_3 - 3 x_4 - x_5)$
$x s.t. \left \{ \begin{array}{l} 2x_2 + x_3 + 4x_4 + 2x_5 \leq 54 \\ 3x_1 + 4x_2 + 5x_3 - x_4 - x_5 \leq 62 \\ x_1, x_2 \geq 0, x_3 \geq 3.32, x_4 \geq 0.678, x_5 \geq 2.57 \end{array} \right.$

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 解析解与数值解

解析解：通过严格推导得到的数学问题的精确解，又称闭式解。
数值解：用数值计算的方法得到的解，通常是问题的近似解。

解析解不存在:$\dfrac{2}{\sqrt{\pi}} \int_0^a e^{-x^2} dx$
数学家的解决方案 – 发明特殊函数:erf(a)
工程技术人员采用近似解。

01.04 计算机数学语言的发展概述

1978 Cleve Moler, New Mexico University, matlab(Matrix Laboratory)。

三个代表: matlab、Mathematica, Maple。

扫描二维码关注公众号，回复： 1785054 查看本文章

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 矩阵乘法通用子程序
$c_{ij} = \sum \limits _{k=1} ^p a_{ik}b_{kj}, i = 1, \cdots, n, j = 1, \cdots, m$

%% matlab 会判断是否满足矩阵A的列数等于矩阵B的行数
C = A * B

01.06 三步求解方法

问题是什么、如何描述、求解。

例1-11a 求函数序列的极限
数学问题：$\lim \limits_{n \to \infty} n \arctan(\dfrac{1}{n(x^2+1) + x}) \tan ^n (\dfrac{\pi}{4} + \dfrac{x}{2n})$
物理意义：
如何描述:

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 线性规划问题求解
$\qquad \min \qquad -2x_1 -x_2 -4x_3 - 3x_4 - x_5$
$x \quad s.t. \left \{ \begin{array}{l}2x_2 + x_3 + 4x_4 + 2x+5 \leq 54 \\ 3x_1 + 4x_2 + 5x_3 - x_4 - x_5 \leq 62 \\ x_1, x_2 \geq 0, x_3 \geq 3.32, x_4 \geq 0.678, x_5 \geq 2.57 \end{array} \right.$
第一步：是什么– 物理含义， 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, ':');