弹簧-质量-阻力系统下图为弹簧-质量-阻力系统的示意图:
根据牛顿力学定理,很容易列出如下的微分方程:
这是一个二阶的微分方程,我们把它改写为如下的两个一阶微分方程组:
下面是一些系统的一些参数,其中外力用f()表示,这样可以用来计算外力随时间变化的情况。
k = 3.0
c = 0.2
m = 0.1
f = lambda t: 10.0
init = 5.0, 0.0
mass_dump_spring()计算在时刻t、状态为status时,和的值。
def mass_dump_spring(status, t):
x, v = status
dx = v
dv = (f(t) - k*x - c*v)/m
return dx, dv
我们可以使用scipy.integrate中提供的odeint()对微分方程进行求解:
import numpy as np
from scipy.integrate import odeint
import pylab as pl
def solve_by_odeint(time, h):
t = np.arange(0, time, h)
result = odeint(mass_dump_spring, init, t)
return t, result[:, 0], result[:, 1]
欧拉方法欧拉方法,命名自它的发明者莱昂哈德·欧拉,是一种一阶数值方法,用以对给定初值的常微分方程(即初值问题)求解。
对于如下的微分方程:
可以通过如下的公式得到时刻的近似值:
下面是实现欧拉方法的程序:
def add(status, dstatus, h): ❶
return [status[i] + dstatus[i]*h for i in xrange(len(status))]
def euler(func, status, time, h):
tlist = np.arange(0, time, h).tolist()
result = []
for t in tlist:
result.append(status)
dstatus = func(status, t)
status = add(status, dstatus, h)
return tlist, np.array(result)
def solve_by_euler(time, h):
t, result = euler(mass_dump_spring, init, time, h)
return t, result[:, 0], result[:, 1]
❶为了支持任意长度的状态矢量,在add()中使用列表推导式计算”status + h * dstatus”的值。此函数在后面实现龙格-库塔法时还会用到。当h足够小时,欧拉方法所得到的值足够精确,但是随着h的增大,误差也会明显增加,下面的程序比较当h为0.01和0.001时的误差:
def euler_plot():
for i, h in enumerate([0.01, 0.001]):
pl.subplot(211 + i)
t, x, v = solve_by_odeint(5, h)
t, x_euler, v_euler = solve_by_euler(5, h)
pl.plot(t, x, label="odeint")
pl.plot(t, x_euler, "r", label="euler")
pl.legend(loc="best")
pl.title("h = %g" % h)
其结果如下图所示:
中点法欧拉法的积累误差与h成正比,为了在较大的h时也能计算出较为精确的解,需要使用更高阶的算法。下面是中点法的计算公式,其积累误差与成正比:
下面是实现中点法的程序,程序中需要调用两次计算微分的函数:
def midpoint(func, status, time, h):
tlist = np.arange(0, time, h).tolist()
result = []
for t in tlist:
result.append(status)
dstatus = func(status, t)
status2 = add(status, dstatus, 0.5*h)
dstatus2 = func(status2, t+0.5*h)
status = add(status, dstatus2, h)
return tlist, np.array(result)
def solve_by_midpoint(time, h):
t, result = midpoint(mass_dump_spring, init, time, h)
return t, result[:, 0], result[:, 1]
经典四阶龙格库塔法龙格库塔法是欧拉方法和中点法的推广,其中最常用的是经典四阶龙格库塔法,通常被称为”RK4”,下面是其计算公式:
其中的等参数使用下面的公式运算,由公式可知需要调用4次计算微分的函数。
下面是实现RK4算法的程序:
def rk4(func, status, time, h):
tlist = np.arange(0, time, h).tolist()
h2 = 0.5*h
result = []
for t in tlist:
result.append(status)
k1 = func(status, t)
k2 = func(add(status, k1, h2), t+h2)
k3 = func(add(status, k2, h2), t+h2)
k4 = func(add(status, k3, h), t+h)
dstatus = [v1+2*v2+2*v3+v4 for (v1,v2,v3,v4) in zip(k1,k2,k3,k4)]
status = add(status, dstatus, h/6)
return tlist, np.array(result)
def solve_by_rk4(time, h):
t, result = rk4(mass_dump_spring, init, time, h)
return t, result[:, 0], result[:, 1]
误差比较为了比较各种算法的积累误差与h之间的关系,我们让h为从0.001到0.1的等比数列,并计算各种算法的运算结果与odeint()的结果之间的误差和。
def error(func1, func2, time, h_list):
ex = []
ev = []
for h in h_list:
_, x1, v1 = func1(time, h)
_, x2, v2 = func2(time, h)
ex.append(np.mean(np.abs(x1-x2)))
ev.append(np.mean(np.abs(v1-v2)))
return ex, ev
def error_plot(func1, func2, title):
h_list = np.logspace(-3, -1, 20)
ex, ev = error(func1, func2, 5.0, h_list)
pl.loglog(h_list, ex, lw=2, label="error x of %s" % title)
pl.loglog(h_list, ev, lw=2, label="error v of %s" % title)
error_plot(solve_by_odeint, solve_by_euler, "euler")
error_plot(solve_by_odeint, solve_by_midpoint, "midpoint")
error_plot(solve_by_odeint, solve_by_rk4, "rk4")
pl.rcParams["legend.fontsize"] = "small"
pl.legend(loc="best")
pl.show()
下图是程序的输出结果,图中X轴和Y轴都是对数坐标,因此在图中误差和h之间成线性关系。而直线的斜率表示算法的阶数。由于odeint()也并非完全精确,其缺省的精度设置为1.49012e-8,因此对于RK4算法,图中部分的误差不再与h成线性关系。
阅读(4261) | 评论(0) | 转发(0) |