一、欧拉法
欧拉法:是一种微分方程的数值计算方法,其基本思想是迭代。所谓迭代,就是逐次替代,最后求出所要求的解,并达到一定的精度。误差可以很容易地计算出来。下面分别介绍欧拉法求解一阶微分方程,二阶微分方程的matlab程序实现。
1.一阶微分方程
{
d
y
(
t
)
d
t
=
f
(
y
(
t
)
,
t
)
y
(
0
)
=
y
0
\begin{cases} {dy(t) \over dt} =f(y(t),t) \\ y(0)=y_0 \end{cases}
{dtdy(t)=f(y(t),t)y(0)=y0
用欧拉法实现的步骤如下:
第一步:在给定的区间[0,T]进行分割,将区间T进行N等分(Δt),分别记为t1,t2,t3,…tn+1;
第二步:根据给定条件,给定一阶导数(导数代表自变量的变化率)和初始值y0,利用talyor展开可以把第n个点t的值y(tn)用f’(tn-1)和tn-1表示出来;
y
(
t
n
+
1
)
=
y
(
t
n
)
+
f
(
y
(
t
)
,
t
)
Δ
t
+
ο
(
Δ
t
2
)
y(t_n+_1) =y(t_n)+f(y(t),t)\Delta t+\omicron(\Delta t^2)
y(tn+1)=y(tn)+f(y(t),t)Δt+ο(Δt2)
第三步:近似,忽略掉高阶无穷小(拉格朗日余项)进行近似,根据一致的y0进行不断迭代;
y
(
t
n
+
1
)
=
y
(
t
n
)
+
f
(
y
(
t
)
,
t
)
Δ
t
y(t_n+_1) =y(t_n)+f(y(t),t)\Delta t
y(tn+1)=y(tn)+f(y(t),t)Δt
Matlab实现如下:
以
{
d
y
(
t
)
d
t
=
e
t
y
(
0
)
=
1
\begin{cases} {dy(t) \over dt} =e^t \\ y(0)=1 \end{cases}
{dtdy(t)=ety(0)=1为例
clear
y0=1;
T=2;
dt=T/N;%分割为1000份(或更多),用于控制精度
t=0:dt:2;
y1=zeros(1,length(t)-1);
y1(1)=y0;
for n=1:length(t)-1
y1(n+1)=y1(n)+exp(t(n))*dt;
end
plot(t,y1,t,exp(t))
可以看出和真实曲线基本完全重合,读者可以利用其他已知函数来验证。
2.二阶微分方程
{
d
2
y
(
t
)
d
t
2
=
f
(
y
(
t
)
,
t
,
y
′
(
t
)
)
y
(
0
)
=
y
0
d
y
(
0
)
d
t
=
v
0
\begin{cases} {d^2y(t) \over dt^2} =f(y(t),t,y'(t)) \\ y(0)=y_0 \\ {dy(0) \over dt} =v_0 \end{cases}
⎩⎪⎨⎪⎧dt2d2y(t)=f(y(t),t,y′(t))y(0)=y0dtdy(0)=v0
具体的原理可以参考‘一阶微分方程’的做法,本系列文章主要用于物理问题的模拟,所以接下来用简谐运动为例来说明,
{
d
2
y
(
t
)
d
t
2
=
−
k
m
y
(
t
)
y
(
0
)
=
y
0
d
y
(
0
)
d
t
=
v
0
\begin{cases} {d^2y(t) \over dt^2} =-\text{\(\frac k m\)} y(t)\\ y(0)=y_0 \\ {dy(0) \over dt} =v_0 \end{cases}
⎩⎪⎨⎪⎧dt2d2y(t)=−mky(t)y(0)=y0dtdy(0)=v0
这个微分方程就有了物理意义,即给定初始位置和速度条件下的简谐运动,可以根据提供的参量通过欧拉法把整个运动定下来,得到任意时刻的速度与位移!
clear
k=1;
m=1;
y0=0;
v0=1;
T=20;
N=100000;%用于控制精度,N越大越接近真实结果
dt=T/N;
t=0:dt:T;
y1=zeros(1,length(t));%产生位移向量
y2=zeros(1,length(t));%产生速度向量
y1(1)=y0;
y2(1)=v0;
for n=1:length(t)-1
y1(n+1)=y1(n)+y2(n)*dt;
y2(n+1)=y2(n)+(-k/m)*y1(n)*dt;
end
plot(t,y1,t,y2)
蓝色线代表y1-t图(位移时间),橙色代表y2-t图(速度时间),读者可以调试k,m,y0,v0得到不同的简谐运动。
相空间
尚若在上边程序中最终执行:
plot(y1,y2)
axis equal%用于使两个坐标轴等刻度
%获得以下图像
通常我们是以x,y为横纵坐标,但也可以用x(位移)和p(动量)分别作为横纵坐标,得到以上图像,即为在相空间中的图像,但是为什么是个正圆呢?
这是因为能量守恒以及k,m的取值决定:
首先在这个例子中没有阻尼,所以能量守恒,即振子的动能和势能合为一个常数,1/2mv^2 + 1/2kx^2=C(常数);
又因为k=m=1,所以上边方程可以化成圆的方程的一般形式所以最终为正圆;
若改变k,m的取值,可想而之最终的结果是个椭圆,读者可以自行尝试。