一阶微分方程基础学习#
Volume Rendering时碰到了微分方程。[THOMAS' CALCULUS -FOURTEENTH EDITION IN SI UNITS]学一下第16章一阶微分方程。1. Solutions, Slope Fields, and Euler’s Method#
一般一阶微分方程的定义和解#
一阶微分方程的形式
\({\Large y' = \frac{dy}{dx} = f(x, y)}\)
例如：
\({\Large y' = x + y}\)
\({\Large y' = x / y}\)
微分方程是要求一个函数，涉及到通用解general solution，总是会带一个常数C。
所以实际有无数个解。
验证某个a(x)是否是一个particular solution。按定义左边求导，右边代入，看是否相等即可。
如果给定一个initial condition(初始条件)，也就是一组(x0,
y0)值，那么可以确定C，得到一个确定的曲线。
而现在要学的是带y的一般形式 \({\Large y'=g(x, y)}\)。
slope fields#
initial condition，并画出此点上的斜率，就得到了slope
fields。见图16.2。那么对于每个微分方程，可以直接画出它的slope fields。 如果再给定一个初始条件。就能顺着slope fields大致找到一条曲线。 这条曲线就是给定初始条件下的微分方程的解！
slope fields可以看出微分方程解的走向，非常有实用价值。
Euler’s Method#
如果我们无法或不需要马上得到一个准确的解，我们可以算出一个近似的数值表(numerical solution)。这叫numerical method。
如果给定微分方程和初始条件
\({\Large y' = \frac{dy}{dx} = f(x, y)}\)
\({\Large y(x_0) = y_0}\)
那么正好可以算线性近似：
\({\Large L(x)=y(x_0) + y'(x_0)(x-x_0)}\)
或
\({\Large L(x)=y(x_0) + f(x_0, y_0)(x-x_0)}\)
线性近似见托马斯微积分第3.9章。
2. First-Order Linear Equations#
一阶线性微分方程的标准形式：
\({\Large \frac{dy}{dx} + p(x)y = q(x)}\)
p和q都是关于x的连续函数。
解Linear Equation#
先把一个微分方程整理成标准形式
\({\Large \frac{dy}{dx} + p(x)y = q(x)}\)
构造一个\(v(x)\)乘到两边
\({\Large v(x)\frac{dy}{dx} + v(x)p(x)y = v(x)q(x)}\)
如果\({\Large v(x)}\)的构造遵循 \({\Large \frac{d}{dx}vy = vy' + vpy}\)
那么可以替代微分方程左边，方程变为
\({\Large \frac{d}{dx}vy = v(x)q(x)}\)
两边都对x积分
\({\Large v(x)y = \int v(x)q(x)dx + C}\)
得到最后的y
\({\LARGE y = \frac{\int v(x)q(x)dx + C}{v(x)} }\)
再看v(x)的构造
从结果
\({\Large \frac{d}{dx}vy = vy' + vpy}\)
反算v
左边按乘法规则
\({\Large v'y + vy' = vy' + vpy}\)
\({\Large v' = vp}\)
\({\Large \frac{dv}{dx} = vp}\)
\({\Large \int\frac{dv}{v} = \int pdx}\)
(这种操作不少见，这里为什么能把求导分离再积分，是有理论依据的，见下一节的说明)
\({\Large ln(v) = \int pdx + C}\)
\({\Large e^{ln(v)} = e^C e^{\int pdx}}\) (两边给e上指数)
得到v的构造
\({\Large v = e^Ce^{\int pdx}}\)
因为v是可以任意构造的，只要符合构造逻辑 \({\Large v' = vp}\) 即可。所以可以选方便的C=0
\({\Large v = e^{\int pdx}}\)
\({\Large v = e^{\int pdx} = e^P}\) (小p积分得到大P，积分后不要C)
\({\Large v' = e^Pp = vp}\) (符合要求)
这些操作对我是重要的，开始没仔细看，后续碰到实际问题很难受，逻辑上总是感觉缺东西。 一般教材里好像也没有详细说明。
一些微积分基础操作#
上一节的疑问，为什么能把求导拆开再积分。
积分中的替换#
见托马斯微积分第5.5章
如果有
\({\Large \int f(g(x))g'(x)dx}\)
可以设\({\Large u=g(x)}\)，把整个积分替换为
\({\Large \int f(u)du}\)
即
\({\Large \int f(u)\frac{du}{dx}dx = \int f(u)du}\)
证明：
设F为f的积分，有
\({\Large \frac{d}{dx}F(g(x)) = f(g(x))g'(x)}\)
两边对x积分，左右换一下，回到问题
\({\Large \int f(g(x))g'(x)dx = \int \frac{d}{dx}F(g(x))dx}\)
\({\Large = F(g(x)) + C}\) (右边F先求导再积出来肯定还是本身)
\({\Large = F(u) + C}\)
\({\Large = \int F'(u)du}\) (先回到关于u的函数，再进积分，就得到了du)
\({\Large = \int f(u)du}\)
反过来操作也成立，即做成展开。
\({\Large \int f(u)du = \int f(u)\frac{du}{dx}dx}\)
一种微分方程的变量分离#
见托马斯微积分第7.4章，或[Calculus Early Transcendental Functions, 7th Edition]6.3章，大多数书都应该会讲。
\({\Large \frac{dy}{dx} = g(x)h(y)}\)
给定这种形式的方程，可以推出一个等式，把x相关和y相关的量各放一边，各自积分结果相等。
\({\Large \int \frac{1}{h(y)}dy = \int g(x)dx}\)
用上一节的替换方法可证明
\({\Large \int \frac{1}{h(y)}dy = \int \frac{1}{h(y)}\frac{dy}{dx}dx}\)
\({\Large = \int \frac{1}{h(y)}g(x)h(y)dx}\) (代入)
\({\Large = \int g(x)dx}\)
3. 应用#
与速度成正比的阻力下的运动#
设物体质量m，距离函数s，速度v，时刻t。
力 = 质量 x 加速度 = \({\Large m\frac{dv}{dt}}\)
阻力与速度成正比
\({\Large m\frac{dv}{dt} = -kv}\)
\({\Large \frac{dv}{dt} = -\frac{k}{m}v}\)
\({\Large \frac{1}{v} \frac{dv}{dt} = -\frac{k}{m}}\) (v!=0，两边可除以v)
\({\Large \int\frac{1}{v} \frac{dv}{dt} dt = \int -\frac{k}{m} dt}\) (两边对t积分)
\({\Large ln(|v|) = -\frac{kt}{m} + C}\)
\({\Large |v| = e^{-\frac{kt}{m} + C}}\)
\({\Large v = \pm e^Ce^{-\frac{kt}{m}}}\)
\({\Large v = Ae^{-\frac{kt}{m}}}\)
其中\({\Large A=\pm e^C}\)
\({\Large v = v_0e^{(-k/m)t}}\)
可以从这个公式看出一些信息。如果质量m很大，那么需要很多时间速度才能降下来。
如果和距离关联
\({\Large \frac{ds}{dt} = v_0e^{(-k/m)t}}\)
两边积分得
\({\Large s = -\frac{v_0m}{k}e^{(-k/m)t} + C}\)
此时可以给定初始条件t=0，s=0。得
\({\Large C = \frac{v_0m}{k}}\)
代入s得
\({\Large s(t) = \frac{v_0m}{k}(1- e^{(-k/m)t})}\)
得到了距离和时间得关系。如果想知道物体最多能滑行多远，t取无穷大。得
\({\Large \lim_{t->\infty}s(t) = \frac{v_0m}{k}(1-0) = \frac{v_0m}{k}}\)
到此如果已知质量和摩擦系数，可以求速度1变化到速度2的时间。可以求某个初始速度下滑行的最大距离。
不准确的人口的指数增长#
\({\Large \frac{dP}{dt} = kP}\)
和上一节的例子相似，假设当前人口变化率为当前总人口乘一个常系数k(相对增长率)。
任意时刻的总人口P为
\({\Large P = P_0e^{kt}}\)
从1980-1989世界人口数据来看，k接近一个常数0.017，而1980年人口为4454m。那么
\({\Large P = P_0e^{kt} = 4454e^{0.017t}}\)
取t=28即算2008年的人口，得7169m。但实际人口为6707m。
为什么不准，因为建模太简单，没考虑其他各种因素。实际上k是会变化的。
4. 几种形式#
对于微分方程的标准形式
\({\Large \frac{dy}{dx} + p(x)y = q(x)}\)
如果q=0且p为常数k，就是上面的摩擦力例子。
\({\Large \frac{dy}{dx} + ky = 0}\)
\({\Large \frac{dy}{dx} = - ky}\)
\({\Large \int \frac{1}{y} dy = \int -k dx}\) (变量分离积分)
\({\Large ln(y) = -kx + C}\)
\({\Large y = e^{- kx + C}}\)
\({\Large y = Ce^{- kx}}\)
如果q=0
\({\Large \frac{dy}{dx} + p(x)y = 0}\)
\({\Large \frac{dy}{dx} = -p(x)y}\)
\({\Large \int \frac{dy}{y} = \int -p(x)dx}\) (变量分离积分)
\({\Large ln(y) = -P(x) + C}\)
\({\Large y = Ce^{-P(x)}}\)
完整形式
\({\LARGE y = \frac{\int v(x)q(x)dx + C}{v(x)} }\)
5.The Fundamental Theorem of Calculus 相关#
微积分基本定理第一部分(见托马斯微积分5.4章)有
\({\Large \frac{d}{dx}\int_{a}^{x} f(t)dt = f(x)}\)
两边对x积分
\({\Large \int_{a}^{x} f(t)dt = \int f(x)dx + C}\)
a可以取任意值。这样可以重写积分。比如微分方程的解
\({\LARGE y = \frac{\int v(x)q(x)dx + C}{v(x)} }\)
可以写成
\({\LARGE y = \frac{\int_{a}^{x} v(t)q(t)dt + C}{v(x)} }\)
这样重写是合法的，可以验证，把v(x)乘回去，两边对x求导，能还原到上一步的
\({\Large \frac{d}{dx}vy = v(x)q(x)}\)
这在某些情况下有用。
此时a仍然可以选定任意值。
中文好像称为”变上限积分”之类。
这里有一个对新手比较重要的思想
\({\Large \int f(x)dx =\int_{a}^{x} f(t)dt + C}\)
就是这样替换之后的积分变量t到底是啥意思，有什么意义？
甚至可以直接不要t，简写成
\({\Large \int f(x)dx =\int_{a}^{x} f + C}\)
另外对于v(x)，也可以转换
\({\Large v(x) = e^{\int pdx} = e^{\int_{a}^{x}p(t')dt' + C} = e^{\int_{a}^{x}p(t')dt'} }\) (同样可以不要C)
同样可以验证
\({\Large v(x) = e^{\int_{a}^{x}p(t')dt'} = e^{P(x) - P(a)} }\)
\({\Large v'(x) = e^{P(x) - P(a)}\frac{d}{dx} (P(x) - P(a)) = e^{P(x) - P(a)} p(x) = vp }\)
不知道是不是正规的写法。\({\Large v(x) = e^{\int_{}^{x}p(t')dt'} = e^{P(x)} }\)
\({\Large v'(x) = e^{P(x)}\frac{d}{dx} P(x) = e^{P(x)} p(x) = vp }\)
然后对于微分方程的解可以做进一步的操作。
因为解的积分变量x换成了t，那么下面的v(x)对于t可看作常数，可以除进去。
\({\LARGE \frac{v(t)}{v(x)} = \frac{e^{\int_{a}^{t}p(t')dt'}}{e^{\int_{a}^{x}p(t')dt'}} = e^{\int_{a}^{t}p(t')dt' - \int_{a}^{x}p(t')dt'} = e^{\int_{x}^{t}p(t')dt'}}\)
那么
\({\LARGE y = \frac{\int_{a}^{x} v(t)q(t)dt + C}{v(x)} }\)
\({\LARGE = \int_{a}^{x} \frac{v(t)}{v(x)}q(t)dt + \frac{C}{v(x)} }\)
\({\LARGE = \int_{a}^{x} e^{\int_{x}^{t}p(t')dt'} q(t)dt + Ce^{-\int_{a}^{x}p(t')dt'} }\)
历尽千辛万苦得到这个结果，它包含一个二重积分。
仍然要记着之前说的，dt和dt’并不重要，积分域和f重要。
\({\LARGE e^{\int_{x}^{t}p(t')dt'}}\)
而此时x，t，p都是已知，求积分，原理上完全没有问题。
做成从t积到x更舒服一些。
\({\LARGE y = \int_{a}^{x} e^{-\int_{t}^{x}p(t')dt'} q(t)dt + Ce^{-\int_{a}^{x}p(t')dt'} }\)
\({\Large y_0=C}\)
(a选0就是为了方便算，以及含义更清晰和通用。你也可以选3选999来算。本质一样，但混乱。)
那么带初始值的解变为
\({\LARGE y(x) = \int_{0}^{x} e^{-\int_{t}^{x}p(t')dt'} q(t)dt + y_0e^{-\int_{0}^{x}p(t')dt'} }\)
其中p,q,y0都是问题的已知数据。
这个实际上是体渲染方程，用来渲染比如烟雾、透光材质等。
例如眼睛在位置0处，看向物体的表面某点，点在x处。眼睛和物体都处于雾霾中。
这个方程就是算这条光路一共对眼睛打过去多少光能，也就是眼睛最后看到的颜色。
右边的y0是特殊的交点，物体的那一点初始能量就是y0，直接算衰减即可。