数学建模算法-day[13]

5.5 函数逼近

如果已知一个较为复杂的连续函数$y(x),x\in [a,b]$，要求选择一个较简单的函数$f(x)$，在一定准则下最接近$y(x)$，就是所谓函数逼近。

与曲线拟合的最小二乘准则对应，函数逼近常用的一种准则是最小平方逼近，即
$$J={\int }_a^b[f(x)-y(x){]}^{2}dx\text{\hspace{1em}(5.15)}$$
达到最小。与曲线拟合一样，选一组函数{rk(x),k=1,2,⋯ ,m}\{r_k(x),k=1,2,\cdots,m\}{rk​(x),k=1,2,⋯,m}构造$f(x)$，即令
$$f(x)=a_1{r}_1(x)+a_2{r}_2(x)+\cdots +a_m{r}_m(x),$$
代入式(5.15),求$a_1,a_2,\cdots ,a_m$使$J$达到极小。利用极值必要条件可得
$$\left[\begin{array}{cccc}\left(r_1,r_1\right)& \left(r_1,r_2\right)& \cdots & \left(r_1,r_m\right)\\ \left(r_2,r_1\right)& \left(r_2,r_2\right)& \cdots & \left(r_2,r_m\right)\\ ⋮& ⋮& \ddots & ⋮\\ \left(r_m,r_1\right)& \left(r_m,r_2\right)& \cdots & \left(r_m,r_m\right)\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_m\end{array}\right]=\left[\begin{array}{c}\left(y,r_1\right)\\ \left(y,r_2\right)\\ ⋮\\ \left(y,r_m\right)\end{array}\right],\text{\hspace{1em}(5.16)}$$

这里 $(g,h)={\int }_a^{b}g(x)h(x)\mathrm{d}x$ 。当方程组（5.16）的系数矩阵非奇异时，有唯一解。
最简单的当然是用多项式逼近函数，即选 $r_1(x)=1,r_2(x)=x,r_3(x)=x^2,\cdots$ 。并且如果能使 ${\int }_a^b{r}_i(x)r_j(x)\mathrm{d}x=0(i\ne j)$ ，方程组（5.16）的系数矩阵将是对角阵，计算大大简化。满足这种性质的多项式称为正交多项式。

勒让得（Legendre）多项式是在 $[-1,1]$ 区间上的正交多项式，它的表达式为

$$P_0(x)=1,P_k(x)=\frac{1}{2^{k}k!}\frac{d^k}{d{x}^k}{\left(x^2-1\right)}^k,k=1,2,\cdots$$

可以证明

$$\begin{array}{c}{\int }_{-1}^1{P}_i(x)P_j(x)\mathrm{d}x=\{\begin{array}{ll}0,& i\ne j,\\ \frac{2}{2i+1},& i=j,\end{array}\\ P_{k+1}(x)=\frac{2k+1}{k+1}x{P}_k(x)-\frac{k}{k+1}{P}_{k-1}(x),k=1,2,\cdots \end{array}$$

常用的正交多项式还有第一类切比雪夫（Chebyshev）多项式

$$T_n(x)=\cos (n\arccos x),x\in [-1,1],n=0,1,2,\cdots$$

和拉盖尔（Laguerre）多项式

$$L_n(x)={\mathrm{e}}^x\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}\left(x^n{\mathrm{e}}^{-x}\right),x\in [0,+\infty ),n=0,1,2,\cdots$$

例 5.21 求 $f(x)=\cos x,x\in \left[-\frac{\pi }{2},\frac{\pi }{2}\right]$ 在 $H=\mathrm{Span}\left\{1,x^2,x^4\right\}$ 中的最佳平方逼近多项式。

matlab程序

% 最佳平方逼近
clc;clear
syms x;
base=[1,x^2,x^4];
A=base'*base;
A=int(A,-pi/2,pi/2);
b=cos(x)*base';
b=int(b,-pi/2,pi/2);
a=double(A\b)

所求的最佳平方逼近多项式为
$$y=0.9996-0.4964x^2+0.0372x^4.$$

其实这章书上还有最后一个比较大的题，就是黄河小浪底跳水调沙问题，额，这个我们在学数值分析的时候做探索型实验做过这个例子，不过当时主播不是负责代码的，主播是负责做ppt的，当时这个例子我们汇报的时候还遭到了老师的疑惑，因为我们第二题弄错了似乎，唉，这个要是主播去做的话那必定是会得到老师的肯定的。那有空再更这个吧，（似乎欠的有点多了…）

6 微分方程

微分方程一般来说都是出现在A题的，且都是非常难的，不仅要列方程，而且还要会解和稳定性分析等等。这类问题对数学和物理能力要求较高，对编程能力则要求不是较高。

微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题，大体上可以按以下步骤进行：

1. 根据实际要求确定要研究的量（自变量、未知函数、必要的参数等）并确定坐标系。
2. 找出这些量所满足的基本规律（物理的、几何的、化学的或生物学的等）。
3. 运用这些规律列出方程和定解条件。

列方程常见的方法有以下几种：

1. 按规律直接列方程。在数学、力学、物理、化学等学科中，许多自然现象所满足的规律已为人们所熟悉，并直接由微分方程所描述，如牛顿第二定律、放射性物质的放射性规律等。我们常利用这些规律对某些实际问题列出微分方程。
2. 微元分析法与任意区域上取积分的方法。自然界中也有许多现象所满足的规律是通过变量的微元之间的关系式来表达的，对于这类问题，我们不能直接列出自变量和未知函数及其变化率之间的关系式，而是通过微元分析法，利用已知的规律建立一些变量（自变量与未知函数）的微元之间的关系式，然后再通过取极限的方法得到微分方程，或等价地通过任意区域上取积分的方法来建立微分方程。
3. 模拟近似法。在生物、经济等学科中，许多现象所满足的规律并不很清楚而且相当复杂，因而需要根据实际资料或大量的实验数据，提出各种假设。在一定的假设下，给出实际现象所满足的规律，然后利用适当的数学方法列出微分方程。

在实际的微分方程建模过程中，也往往是上述方法的综合应用。不论应用哪种方法，通常要根据实际情况，做出一定的假设与简化，并把模型的理论或计算结果与实际情况进行对照验证，以修改模型使之更准确地描述实际问题并进而达到预测预报的目的。

6.1 常微分方程问题的数学模型

例6.1 建立物体冷却过程的数学模型。将物体放置于空气中，在时刻$t=0$时，测得它的温度为$u_0=100^{\circ }C$,20min后测得它的温度为$u_1=60^{\circ }C$。要求建立此物体的温度$u$和时间$t$的关系，并计算经过多长时间此物体的温度将达到$30^{\circ }C$,其中假设空气的温度保持为$20^{\circ }C$。

解 牛顿冷却定理是温度高于周围环境的物体向周围媒质传递热量逐渐冷却时所遵循的规律；当物体表面与周围存在温度差时，单位时间从单位面积散失的热量与温度差成正比，比例系数称为热传递系数，记为$k$.

假设该物体在时刻$t$时的温度为$u=u(t)$，则由牛顿冷却定律，得
$$\frac{du}{dt}=-k(u-20),\text{\hspace{1em}(6.1)}$$
式中:$k>0$。方程(6.1)就是物体冷却过程的数学模型。

可将方程$(6.1)$改写为
$$\frac{d(u-20)}{u-20}=-kdt,$$
两边积分，得
$${\int }_{100}^u\frac{d(u-20)}{u-20}={\int }_0^t-kdt,$$
化简，得
$$u=20+80e^{-kt}.\text{\hspace{1em}(6.2)}$$
将条件$t=20,u=u_1=60$代入式$(6.2)$，得$k=\frac{\ln 2}{20}$，所以此物体的温度$u$和时间$t$的关系为$u=20+80e^{-\frac{\ln 2}{20}t}$.令$30=20+80e^{-\frac{\ln 2}{20}t}$,解之得$t=60$，即60min后物体的温度为$30^{\circ }C$。

例6.2 目标追踪问题。设位于坐标原点的甲舰向位于$x$轴上点$Q_0(1,0)$处的乙舰发射导弹，导弹始终对准乙舰。如果乙舰以最大的速度$v_0(v_0是常数)$沿平行于$y$轴的直线行驶，导弹的速度是$5v_0$，求导弹运行的曲线。又乙舰形式多远时，导弹将它击中？

解 设导弹的轨迹曲线为$y=y(x)$，并设经过时间为$t$,导弹位于点$P(x,y)$，乙舰位于点$Q(1,v_{0}t)$。由于导弹头始终对准乙舰，故此时直线$PQ$就是导弹的轨迹曲线弧$OP$在点$P$处的切线，如图6.1所示

图6.1

则有
$$\frac{dy}{dx}=\frac{v_{0}t-y}{1-x}.\text{\hspace{1em}(6.3)}$$
由于模型中含有参变量$t$，故要求解该模型应增加附加条件。解决这个问题可从我们没有用到的条件$5v_0$下手。

方法1：相同时间的条件下，导弹运动的路程=5*乙舰运动的路程。

故由弧长计算公式可得：
$${\int }_0^x\sqrt{1+(\frac{dy}{dx}{)}^2}dx=5v_{0}t.\text{\hspace{1em}(6.4)}$$
方法2：利用速度分量合成的概念。
$$\sqrt{(\frac{dx}{dt}{)}^2+(\frac{dy}{dt}{)}^2}=5v_0,\text{\hspace{1em}(6.5)}$$
注意到，式(6.5)即式(6.4).

之后要求此问题的定解，还需要给出两个初始条件。事实上，初始时刻轨迹曲线通过坐标顶点，即$x=0,y(0)=0$；此外在该点的切线平行于$x$轴，因此有$y^{\prime }(0)=0$。

对式(6.4)两边求导，得
$$\sqrt{1+(\frac{dy}{dx}{)}^2}=5v_0\frac{dt}{dx}.\text{\hspace{1em}(6.6)}$$

为了消去中间变量$t$，把方程(6.3)改写为
$$(1-x)\frac{dy}{dx}=v_{0}t-y\text{\hspace{1em}(6.7)}$$
然后两边关于$x$求导，整理得
$$(1-x)\frac{d^{2}y}{d{x}^2}=v_0\frac{dt}{dx},\text{\hspace{1em}(6.8)}$$

故联立方程得到：
$$\{\begin{array}{rl}& (1-x)\frac{d^{2}y}{d{x}^2}=v_0\frac{dt}{dx}\\ & \sqrt{1+(\frac{dy}{dx}{)}^2}=5v_0\frac{dt}{dx}.\\ & y(0)=0,y^{\prime }(0)=0.\end{array}$$
消去中间变量$\frac{dt}{dx}$，得关于轨迹曲线的二阶非线性常微分方程：
$$\{\begin{array}{rl}& (1-x)\frac{d^{2}y}{d{x}^2}=\frac{1}{5}\sqrt{1+(\frac{dy}{dx}{)}^2},0<x\le 1\\ & y(0)=0,y^{\prime }(0)=0.\end{array}\text{\hspace{1em}(6.9)}$$
用mathematica软件求得轨迹曲线
$$y(x)=\frac{5}{24}(1-3(1-x{)}^{4/5}+2(1-x{)}^{6/5})\text{\hspace{1em}(6.10)}$$

何处击中乙舰？在式(6.10)中，令$x=1$，得$y=\frac{5}{24}$，，即在$(1,\frac{5}{24})$处击中乙舰。

何时击中乙舰？击中乙舰时，乙舰航行距离$y=\frac{5}{24}$，由$y=v_{0}t$，得$t=\frac{5}{24v_0}$时击中乙舰。

mathematica程序

DSolve[{(1 - x)*y''[x] == (1/5)*Sqrt[1 + (y'[x])^2], y[0] == 0, y'[0] == 0}, y[x], x]

参考文献
[1] 司守奎,孙玺青 数学建模算法与应用第3版[M]

唉，主播晚上看到了一个热点，看的难受的，唉。