Python科学计算：常微分方程计算1

微分方程是一个非常有用的工具，可以说是理工科绕不开的一个坎，刚好，这本书就提到了常微分方程的解法，有些要用到Scipy包，所以，可以先下载一下啊：

pip install scipy

下载好了之后，我们就开始吧，对于一个微分方程来说，我们在初始时刻指定了充分条件以求得确定的解，这种问题被称为初值问题。

今天有可能先不写代码，我们先来看一下初值问题数值积分的基本思想：

我们考虑函数y(t)的单个一阶方程：

那么他的精确解为：

聚焦于初始条件，同时考虑到的非正值，如果，那么精确解就是

对于，解就是0，我们看一看对于取若干非负值时的解：

图1不同lambda

可以看到，随着lambda的绝对值的增大，解会迅速从t=0时的0上升到O(1)，然后以的形式进行衰减，当然，这张图的代码在下面：

import turtle
import matplotlib.pyplot as plt
import math
import numpy as np
y_0=0
t=np.linspace(0.0,1.0,100,dtype=float)
#lambda=-5
x_0=-np.exp(-5*t)+np.exp(-t)
#lambda=-10
x_1=-np.exp(-10*t)+np.exp(-t)
#lambda=-100
x_2=-np.exp(-100*t)+np.exp(-t)
#lambda=-150
x_3=-np.exp(-150*t)+np.exp(-t)
#lambda=-200
x_4=-np.exp(-200*t)+np.exp(-t)
#lambda=-250
x_5=-np.exp(-250*t)+np.exp(-t)
fig=plt.figure()
ax=fig.add_subplot(1,1,1)
ax.plot(t,x_0,"r-",label="lamb=-5")
ax.plot(t,x_1,"b-",label="lamb=-10")
ax.plot(t,x_2,"g-",label="lamb=-100")
ax.plot(t,x_3,"r--",label="lamb=-150")
ax.plot(t,x_4,"b--",label="lamb=-200")
ax.plot(t,x_5,"g--",label="lamb=-250")
ax.legend(loc='best')
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("different lambda")
fig.show()
turtle.done()

其实可以不用写这么多lambda的值的，书上只写了四个，我觉得可以多加两条看看，现在看起来，变化趋势和范围几乎没啥差别。

为了用数值方法来对求解这个方程，引入一个离散网格，选择一些大整数N，并且设置h=1/N,这样，我们就可以定义：

简便起见，选择等距离网格，设，那么在这个网格上面，我们就能把表示成一个序列：,然后令相应的数值近似序列为,

那怎么样才能近似计算呢，最简单的近似防范就是前向欧拉方程：

这种方法，对应于只保留了泰勒级数中的前两项，实际上，截断误差是:

从这里看得出来，如果我们的h选择的够小，那么我们就可以使截断误差尽可能的小(其实，这里我觉得，这种方式你有可能突然接受不了，h是分母啊，怎么h越小截断误差越小，其实这里如果说你看减号后面的还是不甚理解的话，给你个简单的解释，我们常说的微积分，怎么积分才能更加接近结果嘞？微分步长尽可能的小嘛，这样，他就接近了呗)。

截断误差倒是没啥，我们最应该关注的是实际误差：，从上面的式子计算，我们就会得到

对于这种推进关系，我们就可以得到：n>0的时候：

当n时，只有当的时候，才有界，否则的话，就会按照指数形式增长，所以，如果说前向欧拉方程的数值求解结果能够被接受，那么就必须要满足这个条件。

如果说，接近于-1，任何合理的h值都会产生可以接受的结果，但是考虑到，稳定性要求,这在计算快速初始变化时是可以理解的，但是，在其他值域内求解类似缓慢衰减的方程的时候，需要设置非常小的步长，所以数值演化进行的非常缓慢，这种现象被称为问题的刚性；显然，前向欧拉方程在刚性问题上并不令人满意(说实话，没怎么看懂，但结论还是很明确的)。

前向欧拉方具有显式的性质，简单理解，就是已知求未知。

那么怎么样才能解决刚性问题呢？有人就提出了后向欧拉方程的方法：

从这里，可以看到一个很有趣的事情：未知求未知，当然，也可以合并方程中的同类项，其实际误差满足：

其解为：

其稳定性准则为：

很明显，这就满足了我们讨论的对于所有正值h都具有大负值的刚性问题。

后向欧拉方程被称为隐式方案，这些隐式方案主要用来处理刚性问题，