从零开始的微分方程式

2024-02-05 20:55:31 # Math #Math #微分方程式 #数学 #大学

什么是微分方程式？

与一元一次方程可以解出来变量的具体值不同，微分方程式是关于函数的方程，解出来的是一个函数。例如微分方程

$\frac{dy}{dx}=1$ $y=x+C$

其中 $C$ 是积分常数（積分定数）。如果有一个初期条件是 $y(0)=1$ ，即可解得 $C=1,y=x+1$ 。

这篇文章里会从零介绍一階微分方程式、二階微分方程式、偏微分方程式的解法，只需要掌握基础求导和积分的知识。

一/二階微分方程式：方程中最高是求一次或者二次导（微分），例如 $x''+x'+1=3$ 就是一个二階微分方程式。
常微分方程式：方程中只有一个变量，例如变量只有 $x$ 的方程式。
偏微分方程式: 方程中有两个以上的变量，例如包含 $x$ 和 $y$ 的方程组。

为了方便描述， $y'=\frac{dy}{dx}~~or~~\frac{dy}{dt}$

一階微分方程式

前進オイラー法

假设你不会解这个微分方程，只想知道函数 $y(x)$ 在某个点近似的值，可以使用欧拉法近似。

对于每个点 $(x,y)$ 来说，把 $(x,y)$ 带入微分方程式中，可以计算出 $dy/dx$ 的数值，这个值其实就是函数 $y$ 在点 $(x,y)$ 的斜率。

那么可以从设定一个长度 $h$ ，从初始点 $(x_0,y(x_0))$ 开始，计算出横坐标 $x_0+h$ 所对应的 $y$ 的值；如此往复，直到计算出目标值 $y(x)$ 。

$y(x_0+h)=y(x_0)+h\times \frac{dy}{dx} |_{x=x_0,y=y_0}$

显然，$h$的值越小，需要的计算次数越多，得到的结果越精确。

例：对于微分方程式

$\frac{dy}{dx}=-\frac{y}{x}+3x$

初期条件 $y(1)=2$ ，令 $h=0.1$ ，计算 $y(1.3)$ 的值。

解：可以画出下图的表格，方便计算。

代码实现

#include <bits/stdc++.h>
using namespace std;
const int maxn = 1e5 + 10;
double x[maxn], y[maxn], A[maxn], h = 0.1;
double F(double x, double y){
    return - (y / x) + 3 * x;    
}
int main(){
    x[0] = 1.0, y[0] = 4.0, A[0] = 0;
    printf("y(?)=?\n");
    scanf("%lf%lf",&x[0],&y[0]);
    printf("h=?\n");
    scanf("%lf",&h);
    printf("ask for y(?)\n");
    double End;
    scanf("%lf", &End);
    A[0] = F(x[0], y[0]);
    int times = int((End - x[0]) / h);
    for (int i = 1; i <= times; ++i) {
        x[i] = x[i-1] + h;
        y[i] = y[i - 1] + A[i - 1] * h;
        A[i] = F(x[i], y[i]); 
    }
    printf("-------------------------------------\n");
    for (int i = 0; i <= times; ++i)
        printf("%d %lf %lf %lf %lf\n", i, x[i], y[i], A[i], h * A[i]);
    printf("y(%lf)=%lf", End, y[times]);
}

変数分離法

对于方程

$\frac{dy}{dx}=f(x)g(y)\\$

我们可以把y和x有关的项分离，放到等号的两边，然后积分。

假设 $F(x)'=f(x)$

$\begin{align} \frac{1}{g(y)}{dy}&=f(x)dx\\ \int \frac{1}{g(y)}{dy}&=\int f(x)dx\\ \ln(g(y))&=F(x)+C\\ g(y)&=e^{F(x)}+e^C\\ g(y)&=C'e^{F(x)} ~~~~(C'=e^C) \end{align}$

例：

$\begin{align} \frac{dy}{dx}+3y&=0\\ \end{align}$

解：

$\begin{align} \frac{dy}{dx}&=-3y\\ \frac{1}{y}{dy}&=-3dx\\ \ln(y)&=-3x+C\\ y&=C'e^{-3x} \end{align}$

同次形

如果 $x$ 和 $y$ 的次数一致，可以用 $u$ 替换 $\frac{y}{x}$ ，把方程化简以后可以使用変数分離法求解。

例：

$\begin{align} (xy+y^2+x^2)dx-x^2dy&=0\\ x^2dy&=(xy+y^2+x^2)dx\\ \frac{dy}{dx}&=(\frac{y}{x}+\frac{y^2}{x^2}+1)\\ \end{align}$

这里引入一个新变量 $u$

$\begin{align} let~~u&=\frac{y}{x}\\ \because y&=ux\\ \therefore y'&=u'x+u \end{align}$

带入方程中

$\begin{align} \frac{du}{dx}x+u&=u+u^2+1\\ \frac{du}{dx}x&=u^2+1\\ \frac{du}{dx}&=\frac{u^2+1}{x}\\ \frac{du}{u^2+1}&=\frac{dx}{x}\\ \tan^{-1}(u)&=\ln(x)+C\\ u&=\tan(ln(x)+C)\\ y&=ux= \tan(ln(x)+C)x \end{align}$

積分因子法

对于以下形状的微分方程式

$y'+P(x)y=Q(x)\\$

这里，积分因子u(x)定义为 $e^{\int P(x)dx}$

$\begin{align} u'(x)&=(e^{\int P(x)dx})'\\&=e^{\int P(x)dx}(\int P(x)dx)'\\&=e^{\int P(x)dx}P(x)\\&=u(x)P(x) \end{align}$

如果把方程两边乘上积分因子，发现方程式左边其实变成了 $uy$ 在微分以后的形式

$\begin{align} uy'+uP(x)y&=uQ(x)\\ (uy)'&=uQ(x)\\ \int (uy)'dx&=\int uQ(x)dx\\ uy&=\int uQ(x)dx+C\\ y&=u^{-1}\int uQ(x)dx+C\\ \end{align}$

例，在变数分离法的例题上用积分因子法，也能得到同样的结果

$\begin{align} \frac{dy}{dx}+3y&=0\\ \end{align}$

解：

$\begin{align} \frac{dy}{dx}+3y&=0\\ e^3\frac{dy}{dx}+3e^3y&=0\\ (e^3y)'&=0\\ e^3y&=C\\ y&=Ce^{-3x} \end{align}$

参考文献

真貝寿明. (2010). 徹底攻略常微分方程式.
微分方程式の数値解法. 東京大学工学部精密工学科プログラミング応用 I・ II. http://www.den.t.u-tokyo.ac.jp/ad_prog/ode/
予備校のノリで学ぶ「大学の数学・物理」. (2020).【大学数学】微分方程式入門 [Vedio]. Youtube. https://www.youtube.com/watch?v=po97dnBfoco
NHCLDX. (2019). 第七章——微分方程. https://zhuanlan.zhihu.com/p/99271908

2024-02-05 20:55:31 # Math #Math #微分方程式 #数学 #大学