欧拉法求解微分方程
Posted conmajia
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了欧拉法求解微分方程相关的知识,希望对你有一定的参考价值。
by Conmajia
2014
工程中有时候需要解微分方程,比如这种:
[ y'=y-cfrac{2x}{y} ]
(y(0)=1). 两边积分,得到:
[ egin{align*} int_{x_n}^{x_{n+1}}y'mathrm{d} x&=int_{x_n}^{x_{n+1}}left(y-cfrac{2x}{y} ight)mathrm{d} x y_{n+1}-y_n&=int_{x_n}^{x_{n+1}}left(y-cfrac{2x}{y} ight)mathrm{d}x y&=sqrt{1+2x} end{align*} ]
工程上只想要数值解,一般采用差分近似代替微分. 最简单最朴素的办法是用欧拉公式:
[ ag{1} y_{n+1}=y_n+Delta x f(x_n,y_n)quad n=0,1,2,cdots ]
推导很简单,对(x_n),有:
[ y'_n=f(x_n,y_n). ]
对(y'(x_n))有:
[ ag{2} y'_napprox cfrac{y_{n+1}-y_n}{Delta x}. ]
式(2)左边叫微商(不是喜提迪丽热巴的那个微商),右侧叫差商(不是淘宝伪劣产品卖家),(Delta x=left|x_{n+1}-x_n ight|).
代回式(2),有:
[ egin{align*} y'_n=f(x_n,y_n) Rightarrowcfrac{y_{n+1}-y_n}{Delta x} &approx f(x_n,y_n) y_{n+1}-y_n &approx Delta x f(x_n,y_n) y_{n+1} &approx y(x_n)+Delta x f(x_n,y_n). end{align*} ]
式(1)的欧拉公式成了:
[ ag{3} y_{n+1}=y_n+Delta xleft(y_n-cfrac{2x_n}{y_n} ight). ]
假设用myFn
函数表示余项(y_n-dfrac{2x_n}{y_n}):
function myFn(x, y) {
return y - 2 * x / y;
}
Euler可以这样实现:
function calculate(x0, y0, delta, xn) {
var yn;
while(x0 < xn) {
yn = y0 + delta * myFn(x0, y0);
y0 = yn;
x0 = x0 + delta;
}
return yn;
}
现在可以开始试验了。
理论上:
[ y=sqrt{1+2x}Rightarrow y(1)=sqrt{3}approx 1.7321 ]
(sqrt{3})是方程的真值,程序的目标是通过计算,得到尽量接近真值的结果.
取(x_0=0),(y_0=1),(Delta x=0.1),程序计算结果为:
ans = 1.784771
误差3.04%. 减小(Delta x),比如(Delta x=0.0001),计算结果为:
ans = 1.732112
误差0.0007%,十分接近真值了.
虽然这个方法可以求到比较精确的解,但是(Delta x)太小的话,while
会执行很多次,效率低下.
引入定积分的梯形公式:
[ int_{x_n}^{x_{n+1}}f(x,y)mathrm{d}xapproxcfrac{Delta x}{2}left[f(x_n,y_n)+f(x_{n+1},y_{n+1}) ight] ]
式(3)可以写成:
[ ag{4} y_{n+1}approx y_n+cfrac{Delta x}{2}left[f(x_n,y_n)+f(x_{n+1},y_{n+1}) ight]. ]
式(3)和式(4)的特点是,前者速度快,精度低;后者速度慢,精度高. 所以对于粗算,可以使用:
[ y'_{n+1}=y_n+Delta x f(x_n,y_n) ]
精算,可以用:
[ y_{n+1}=y_n+cfrac{Delta x}{2}left[f(x_n,y_n)+f(x_{n+1},y'_{n+1}) ight]. ]
结合起来,先通过粗算得到(y'_{n+1})的近似值,再进行精算,得到高精度的最终解. 这样的话,既保证了计算结果的准确度,又没有消耗太多的计算资源,保证了计算效率. 下面是改进后的计算程序:
function calculate(x0, y0, delta, xn) {
var yp, yc;
while(x0 < xn) {
yp = y0 + delta * myFn(x0, y0);
yc = y0 + delta * myFn(x0 + delta, yp);
y0 = 1 / 2 * (yp + yc);
x0 = x0 + delta;
}
return y0;
}
取(x_0=0),(y_0=1),(Delta x=0.1),计算结果为:
ans = 1.737867
误差0.33%.
和最早版本(误差3.04%)相比,在(Delta x)相同——意味着循环次数相近——的情况下,精度提高接近10倍.
具体问题具体分析,手工求微分方程基本是这么个思路. 当然,对于《数值分析》和《工程数学》这些课程来说,我上面写的东西不过是小儿科了. 更多的时候,做个伸手党其实很不错的,大把的现成玩意儿可以用,我干嘛还要费那劲呢?
The End. (Box)
以上是关于欧拉法求解微分方程的主要内容,如果未能解决你的问题,请参考以下文章
MATLAB-1个实例-欧拉法改进欧拉法ode45求解微分方程
数值分析:用改进欧拉法解微分方程初值问题(vf编程) 100