数值方法求解微分方程

一、欧拉法

欧拉法：是一种微分方程的数值计算方法，其基本思想是迭代。所谓迭代，就是逐次替代，最后求出所要求的解，并达到一定的精度。误差可以很容易地计算出来。下面分别介绍欧拉法求解一阶微分方程，二阶微分方程的matlab程序实现。

1.一阶微分方程

$$\{\begin{array}{l}\frac{dy(t)}{dt}=f(y(t),t)\\ y(0)=y_0\end{array}$$
用欧拉法实现的步骤如下：

第一步：在给定的区间[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)\Delta t+o(\Delta t^2)$$

第三步：近似，忽略掉高阶无穷小（拉格朗日余项）进行近似，根据一致的y0进行不断迭代；
$$y(t_n{+}_1)=y(t_n)+f(y(t),t)\Delta t$$
Matlab实现如下：
以$$\{\begin{array}{l}\frac{dy(t)}{dt}=e^t\\ y(0)=1\end{array}$$为例

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.二阶微分方程

$$\{\begin{array}{l}\frac{d^{2}y(t)}{d{t}^2}=f(y(t),t,y^{\prime }(t))\\ y(0)=y_0\\ \frac{dy(0)}{dt}=v_0\end{array}$$
具体的原理可以参考‘一阶微分方程’的做法，本系列文章主要用于物理问题的模拟，所以接下来用简谐运动为例来说明，
$$\{\begin{array}{l}\frac{d^{2}y(t)}{d{t}^2}=-\text{mk}y(t)\\ y(0)=y_0\\ \frac{dy(0)}{dt}=v_0\end{array}$$
这个微分方程就有了物理意义，即给定初始位置和速度条件下的简谐运动，可以根据提供的参量通过欧拉法把整个运动定下来，得到任意时刻的速度与位移！

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的取值，可想而之最终的结果是个椭圆，读者可以自行尝试。

总结：未完待续！