作者自我介绍：大爽歌, b站小UP主 ，直播编程+红警三 ，python1对1辅导老师 。
本博客为一对一辅导学生python代码的教案, 获得学生允许公开。
具体辅导内容为《Computational Problems for Physics With Guided Solutions Using Python》，第一章（Chapter 1 Computational Basics for Physics）的 Code Listings
由于代码较多，第一章内容分为三部分
A(1.1 - 1.6)
B(1.7 - 1.11)
C(1.12 - 1.17)
1 - 1.14 rk4Algor.py
1.12的 rk4Call.py需要导入1.14的rk4Algor.py，
所以先介绍下rk4Algor.py。
1 代码展示
# rk4Algor.py
def rk4Algor(t, h, N, y, f):
k1 = [0]*(N); k2 = [0]*(N); k3 = [0]*(N); k4 = [0]*(N)
fvector = [0]*(N); ydumb = [0]*(N)
fvector = f(t, y) # Returns RHS's
for i in range(0, N): k1[i] = h*fvector[i]
for i in range(0, N): ydumb[i] = y[i] + k1[i]/2.
fvector = f(t+h/2., ydumb)
for i in range(0, N): k2[i] = h*fvector[i]
for i in range(0, N): ydumb[i] = y[i] + k2[i]/2.
fvector = f(t+h/2., ydumb)
for i in range(0, N): k3[i] = h*fvector[i]
for i in range(0, N): ydumb[i] = y[i] + k3[i]
fvector = f(t+h, ydumb)
for i in range(0, N): k4[i] = h*fvector[i]
for i in range(0, N): y[i] = y[i] + (k1[i]+2.*(k2[i]+k3[i])+k4[i])/6.
return y
2 公式原理
doc: Runge-Kutta 4th Order Method to Solve Differential Equation
公式具体说明见链接，这里只分析其计算方法。
使用代码中对应的变量，计算公式写法如下：
k
1
=
h
f
(
t
,
y
n
)
k_1 = hf(t, y_n)
k1=hf(t,yn)
k
2
=
h
f
(
t
+
h
2
,
y
n
+
k
1
2
)
k_2 = hf(t+\frac{h}{2}, y_n+\frac{k_1}{2})
k2=hf(t+2h,yn+2k1)
k
3
=
h
f
(
t
+
h
2
,
y
n
+
k
2
2
)
k_3 = hf(t+\frac{h}{2}, y_n+\frac{k_2}{2})
k3=hf(t+2h,yn+2k2)
k
4
=
h
f
(
t
+
h
2
,
y
n
+
k
3
2
)
k_4 = hf(t+\frac{h}{2}, y_n+\frac{k_3}{2})
k4=hf(t+2h,yn+2k3)
y
n
+
1
=
y
n
+
k
1
6
+
k
2
3
+
k
3
3
+
k
4
6
+
O
(
h
5
)
y_{n+1} = y_n + \frac{k_1}{6} + \frac{k_2}{3} + \frac{k_3}{3} + \frac{k_4}{6} + O(h^5)
yn+1=yn+6k1+3k2+3k3+6k4+O(h5)
其
中
O
(
h
5
)
代
表
偏
差
量
，
使
用
公
式
时
不
计
算
，
好
像
也
算
不
了
其中O(h^5)代表偏差量，使用公式时不计算，好像也算不了
其中O(h5)代表偏差量，使用公式时不计算，好像也算不了
3 逐行解析
解析如下图所示
在每一个阶段，ydumb起的是阶段内部存储暂时值的作用，
4 运行示例
以后面将要调用rk4Algor方法的rk4Call.py中的一次调用的参数为例
- t: 0.0
- h: 0.1
- y: array[ 3. -5.]
- f: f ( t , y ) = [ y 1 , − 100 y 0 − 2 y 1 + 10 s i n ∗ ( 3 t ) ] f(t, y) = [y_1, -100y_0 - 2y_1 + 10sin*(3t)] f(t,y)=[y1,−100y0−2y1+10sin∗(3t)]
其运行后的结果（返回的y）如下
y: array[ 1.34656554 -24.67337214]
其中间过程为
k1: [-0.5, -29.0]
k2: [-1.9500000000000002, -23.450561867526403]
k3: [-1.6725280933763202, -16.755505680773762]
k4: [-2.175550568077376, -8.628097723420707]
2 - 1.12 rk4Call.py
1 代码展示
# rk4Call.py
from vpython import *
from rk4Algor import rk4Algor # 1.14 rk4Algor.py
import numpy as np
Tstart = 0.; Tend = 10.; Nsteps = 100 # Initialization
y = np.zeros((2), float)
graph1 = canvas(x=0, y=0, width=400, height=400, title='RK4',
xtitle='t', ytitle='y[0] = Position versus Time', xmin=0,
xmax=10, ymin=-2, ymax=3)
funct1 = gcurve(color=color.yellow)
graph2 = canvas(x=400, y=0, width=400, height=400, title='RK4',
xtitle='t', ytitle='y[1] = Velocity versus Time',
xmin=0, xmax=10, ymin=-25, ymax=18)
funct2 = gcurve(color=color.red)
y[0] = 3.; y[1] = -5. # Initial position & velocity
t = Tstart; h = (Tend - Tstart) / Nsteps;
def f(t, y): # Force (RHS) function
fvector = np.zeros((2), float)
fvector[0] = y[1]
fvector[1] = -100. * y[0] - 2. * y[1] + 10. * sin(3. * t)
return fvector
if ((t + h) > Tend): h = Tend - t # Last step
y = rk4Algor(t, h, 2, y, f)
t = t + h
rate(30)
funct1.plot(pos=(t, y[0]))
funct2.plot(pos=(t, y[1]))
2 绘制效果
这个绘制效果是一个从左到右绘制图表的动画，。
3 np.zeros
基础使用方法为
np.zeros(shape, dtype=float)
doc: numpy.zeros
Return a new array of given shape and type, filled with zeros.
返回给定维度和类型的新数组，数组中填充0。
常用于创建初始的数组。
- shape: int or tuple of ints.
新数组的维度。
一个整数代表一维（整数值代表元素数量），
元组代表多维（元组的一个元素代表一个维，其值代表这个维的元素数量）， - dtype: data-type, optional.
数组所需的数据类型。默认为浮点数np.float64
使用示例(来自官方文档)
>>> np.zeros(5)
array([ 0., 0., 0., 0., 0.])
>>> np.zeros((5,), dtype=int)
array([0, 0, 0, 0, 0])
>>> np.zeros((2, 1))
array([[ 0.],
[ 0.]])
>>> s = (2,2)
>>> np.zeros(s)
array([[ 0., 0.],
[ 0., 0.]])
4 逐行解析
Tstart = 0.; Tend = 10.; Nsteps = 100 # Initialization
- Tstart: 动画开始时间
- Tend: 动画结束时间
- Nsteps: 动画分多少步（段），对应动画有多少关键帧。
y = np.zeros((2), float)
创建一个一维二元数组，元素都为0：array([0., 0.])
graph1 = graph(x=0,y=0,width=400,height=400,title='RK4',
xtitle='t', ytitle='y[0] = Position versus Time',xmin=0,
xmax=10,ymin=-2,ymax=3)
funct1 = gcurve(color = color.yellow)
graph2 = graph(x=400,y=0,width=400,height=400,title='RK4',
xtitle='t', ytitle='y[1] = Velocity versus Time',
xmin=0,xmax=10, ymin=-25,ymax=18)
funct2 = gcurve(color = color.red)
- graph1, graph2: 两个图,用于分别展示两个曲线。
xmin, xmax: x轴的范围，
ymin, ymax: y轴的范围。 - funct1, funct2: 定义俩个曲线，一个黄色，一个红色
y[0] = 3.; y[1] = -5. # Initial position & velocity
t = Tstart; h = (Tend-Tstart)/Nsteps;
初始化一些值。
h: 一帧的时间长度。计算出来结果为0.1
def f(t, y): # Force (RHS) function
fvector = np.zeros((2), float)
fvector[0] = y[1]
fvector[1] = -100.*y[0]-2.*y[1] + 10.*sin(3.*t)
return fvector
- 定义一个一维二元数组
fvector - 第一个元素为 y 1 y_1 y1
- 第二个元素为 − 100 y 0 − 2 y 1 + 10 s i n ∗ ( 3 t ) -100y_0 - 2y_1 + 10sin*(3t) −100y0−2y1+10sin∗(3t)
while (t < Tend): # Time loop
if ((t + h) > Tend): h = Tend - t
进行时间循环，循环持续到终止时间
如果当前时间+下一帧时间长度大于终止时间，修改下一帧时间长度。
y = rk4Algor(t, h, 2, y, f)
t = t + h
rate(30)
funct1.plot(pos = (t, y[0]) )
funct2.plot(pos = (t, y[1]) )
调用rk4Algor方法，计算出y。
进入下一帧。
设置循环间隔为1/30秒，每秒t增长30*0.1=3，相当于三倍速展示。
- 曲线1添加点 ( t , y 0 ) (t, y_0) (t,y0)
- 曲线2添加点 ( t , y 0 ) (t, y_0) (t,y0)
3 - 1.13 rk4Duffing.py
1 代码展示
import numpy as np, matplotlib.pylab as plt
from math import *; from rk4Algor import rk4Algor
tt = np.zeros((2000),float); yy = np.zeros((2000),float);
vy = np.zeros((2000),float); y = [0]*(2); rhs = [0]*(2)
a = 0.5; b = -0.5; g = 0.02; F = 0.0008; w = 1.
h = 0.1 # Time step
i = 0; y[0] = 0.09; y[1] = 0 # Initial x, velocity
def f(t,y): # RHS function
rhs[0] = y[1]
rhs[1] = -2*g*y[1] - a*y[0] - b*y[0]**3 + F*cos(w*t)
return rhs
for t in np.arange(0,200,h): # Time Loop
tt[i] = t
y = rk4Algor(t,h,2,y,f) # Call rk4
yy[i] = y[0] # x(t)
vy[i] = y[1] # v(t)
i = i+1
fig, axes = plt.subplots(nrows=1, ncols=2,figsize=(12,5) )
axes[0].plot(tt[1000:],yy[1000:]) # 1000 to avoid transients
axes[0].grid() # x(t)
axes[0].set_title('Duffing Oscillator x(t)')
axes[0].set_xlabel('t')
axes[0].set_ylabel('x(t)')
axes[1].plot(yy[1000:],vy[1000:])
axes[1].grid()
axes[1].set_title('Phase Space Orbits for Duffing Oscillator')
axes[1].set_xlabel('x(t)')
axes[1].set_ylabel('v(t)')
plt.show()
2 绘制效果
3 逐行分析
tt = np.zeros((2000),float); yy = np.zeros((2000),float);
vy = np.zeros((2000),float); y = [0]*(2); rhs = [0]*(2)
a = 0.5; b = -0.5; g = 0.02; F = 0.0008; w = 1.
h = 0.1 # Time step
i = 0; y[0] = 0.09; y[1] = 0 # Initial x, velocity
初始化
- tt, yy, vy: 一维2000元的0.0填充的数组。
- y, rhs: 有两个0组成的列表。
- a, b, g, F, w: rhs函数中用到的参数。
- h: 最小的时间间距，两个相邻状态之间的时间间隔。
- i: 计数，序号。
- y: 设置为[0.09, 0]
def f(t,y): # RHS function
rhs[0] = y[1]
rhs[1] = -2*g*y[1] - a*y[0] - b*y[0]**3 + F*cos(w*t)
return rhs
新建函数
f
(
t
,
y
)
=
[
y
1
,
−
2
g
y
1
−
a
y
0
−
b
y
0
3
+
F
c
o
s
(
w
t
)
]
f(t, y) = [y_1, -2gy_1 - ay_0-by_0^3 + Fcos(wt)]
f(t,y)=[y1,−2gy1−ay0−by03+Fcos(wt)]
for t in np.arange(0,200,h): # Time Loop
tt[i] = t
y = rk4Algor(t,h,2,y,f) # Call rk4
yy[i] = y[0] # x(t)
vy[i] = y[1] # v(t)
i = i+1
- 在0到200之间，隔h(0.1)取值，最终能取到2000个值。
- tt数组记录每一个时间值
- 调用
rk4Algor计算出y - yy, vy: 分别记录y的两个值。
- 计数器i加一
fig, axes = plt.subplots(nrows=1, ncols=2,figsize=(12,5))
划分窗口为一行两列个图
- figsize: 指定图形尺寸(宽、高)，单位为英寸。
axes[0].plot(tt[1000:],yy[1000:]) # 1000 to avoid transients
axes[0].grid() # x(t)
axes[0].set_title('Duffing Oscillator x(t)')
axes[0].set_xlabel('t')
axes[0].set_ylabel('x(t)')
axes[1].plot(yy[1000:],vy[1000:])
axes[1].grid()
axes[1].set_title('Phase Space Orbits for Duffing Oscillator')
axes[1].set_xlabel('x(t)')
axes[1].set_ylabel('v(t)')
plt.show()
在子图1中绘制tt[1000:]和yy[1000:]曲线。
在子图2中绘制yy[1000:]和vy[1000:]曲线。
4 1.15 TrapMethods.py
1 代码与运行结果
from numpy import *
def func(x):
return 5*(sin(8*x))**2*exp(-x*x)-13*cos(3*x)
def trapezoid(A,B,N):
h = (B - A)/(N - 1) # step size
sum = (func(A)+func(B))/2 # (1st + last)/2
for i in range(1, N-1):
sum += func(A+i*h)
return h*sum
A = 0.5
B = 2.3
N = 1200
print(trapezoid(A,B,N-1))
运行输出如下
2.998428176678047
2 逐行解析
def func(x):
return 5*(sin(8*x))**2*exp(-x*x)-13*cos(3*x)
定义函数
f
(
x
)
=
5
s
i
n
2
(
8
x
)
e
−
x
2
−
13
c
o
s
(
3
x
)
f(x) = 5sin^2(8x)e^{-x^2} - 13cos(3x)
f(x)=5sin2(8x)e−x2−13cos(3x)
def trapezoid(A,B,N):
h = (B - A)/(N - 1) # step size
sum = (func(A)+func(B))/2 # (1st + last)/2
for i in range(1, N-1):
sum += func(A+i*h)
return h*sum
定义函数，内部执行的代码对应公式如下
h
=
B
−
A
N
−
1
s
t
r
a
p
e
z
o
i
d
(
A
,
B
,
N
)
=
h
(
f
(
A
)
+
f
(
B
)
2
+
∑
i
=
1
N
−
1
f
(
A
+
i
h
)
)
\begin{aligned} h &= \frac{B - A}{N - 1} \\s trapezoid(A,B,N) &= h(\frac{f(A)+ f(B)}{2} + \sum_{i=1}^{N-1} f(A+ih)) \end{aligned}
hstrapezoid(A,B,N)=N−1B−A=h(2f(A)+f(B)+i=1∑N−1f(A+ih))
该公式为梯形积分公式，参考文档：trapezoidal-rule
5 1.17 GaussPoints.py
1.16的 IntegGaussCall.py需要导入1.17的GaussPoints.py，
所以先介绍下GaussPoints.py。
1 代码展示
import math
from numpy import *
def GaussPoints(Npts, a, b, x, w, eps):
m = 0; i = 0; j = 0; t = 0.; t1 = 0.; pp = 0.
p1 = 0.; p2 = 0.; p3 = 0.
m = int((Npts+1)/2)
for i in range(1, m+1):
t = math.cos(math.pi*(float(i)-0.25)/(float(Npts)+0.5))
t1 = 1
while((abs(t-t1)) >= eps):
p1 = 1. ; p2 = 0.
for j in range(1, Npts + 1):
p3 = p2; p2 = p1
p1 = ((2.*float(j)-1)*t*p2 - (float(j)-1.)*p3)/(float(j))
pp = Npts*(t*p1 - p2)/(t*t - 1.)
t1 = t
t = t1 - p1/pp
x[i-1] = -t
x[Npts-i] = t
w[i-1] = 2./( (1.-t*t)*pp*pp)
w[Npts-i] = w[i-1]
for j in range(0, Npts): # Scale [-1,+1] to [a,b]
x[j] = x[j]*(b-a)/2. + (b+a)/2.
w[j] = w[j]*(b-a)/2.
2 公式原理
好像是高斯-勒让得求积
Gauss-Legendre Quadrature
公式有点复杂，没搞太明白
相关文章链接：
- Legendre-Gauss Quadrature
- Gauss-Legendre_Quadrature
- wiki/Gauss–Legendre_quadrature
- wiki/Legendre_polynomials
3 逐行解析
def GaussPoints(Npts, a, b, x, w, eps):
定义函数名，函数参数为
- Npts: 后面用n表示
- a, b: 区间
- x, w: 存放值的数组
- eps: 精度要求。
m = 0; i = 0; j = 0; t = 0.; t1 = 0.; pp = 0.
p1 = 0.; p2 = 0.; p3 = 0.
声明了一堆变量，但其实都是多余了声明，除了pp。
m = int((Npts+1)/2)
for i in range(1, m+1):
m 是N的一半取整(int是向下取整)
因为x,w首尾都是对称的，所以只需循环N的一半。
t = math.cos(math.pi*(float(i)-0.25)/(float(Npts)+0.5))
t1 = 1
对应公式
t
=
c
o
s
(
π
i
−
1
4
n
+
1
2
)
t = cos(\pi \frac{i - \frac{1}{4}}{n + \frac{1}{2}})
t=cos(πn+21i−41)
来自直播间老哥的补充：猜测是和勒让得(Legendre)多项式是有关。
while((abs(t-t1)) >= eps):
不断循环，直到t和t1的差值小于eps
p1 = 1. ; p2 = 0.
for j in range(1, Npts + 1):
p3 = p2; p2 = p1
p1 = ((2.*float(j)-1)*t*p2 - (float(j)-1.)*p3)/(float(j))
先看for循环，会发现
- p2 取上一轮的p1的值
- p3 取上一轮的p2的值，也就是上上一轮的p1的值。
然后用上一轮的p1和上上一轮的p1 不断的计算新的p1，计算N遍。
把这个转化成数列，则如下
数
列
{
p
n
}
,
p
0
=
1
,
p
−
1
=
0
,
p
j
=
(
2
j
−
1
)
t
p
j
−
1
−
(
j
−
1
)
p
(
j
−
2
)
j
\begin{aligned} 数列\{p_n\}, p_0 = 1, p_{-1} = 0, p_j = \frac{(2j-1)tp_{j-1} - (j-1)p_(j-2)}{j} \end{aligned}
数列{pn},p0=1,p−1=0,pj=j(2j−1)tpj−1−(j−1)p(j−2)
执行完for循环(共N轮)后，最后的p1即pN,p2即p(N-1)
来自直播间老哥的补充：猜测是勒让得(Legendre)的递推式。
pp = Npts*(t*p1 - p2)/(t*t - 1.)
t1 = t
t = t1 - p1/pp
pp是 N t p N − p N − 1 t 2 − 1 N\frac{tp_N - p_{N-1}}{t^2 - 1} Nt2−1tpN−pN−1
来自直播间老哥的补充：猜测是勒让得(Legendre)的求导式。
t1 取t的值
t自减pN/pp
x[i-1] = -t
x[Npts-i] = t
w[i-1] = 2./( (1.-t*t)*pp*pp)
w[Npts-i] = w[i-1]
x数组的第i个和倒数第i个值，分别取值-t和t
w数组的第i个和倒数第i个值，都取值
2
(
1
−
t
2
)
p
p
2
\frac{2}{(1-t^2)pp^2}
(1−t2)pp22和
2
(
1
−
t
2
)
p
p
2
\frac{2}{(1-t^2)pp^2}
(1−t2)pp22
for j in range(0, Npts): # Scale [-1,+1] to [a,b]
x[j] = x[j]*(b-a)/2. + (b+a)/2.
w[j] = w[j]*(b-a)/2.
缩放x, w两个数组的值， 从[-1,+1] 缩放为 [a,b]
由于x, w都是数组，其实可以不用for循环，直接对数组进行公示计算，即可改为
x = x * (b - a) / 2. + (b + a) / 2.
w = w*(b-a)/2.
6 1.16 IntegGaussCall.py
1 代码与运行结果
from numpy import *; from GaussPoints import GaussPoints
Npts = 10; Ans = 0; a = 0.; b = 1.; eps = 3.E-14
w = zeros(2001, float); x = zeros(2001, float) # Arrays
def f(x): return exp(x) # Integrand
GaussPoints(Npts, a, b, x, w, eps) # eps: precison of pts
for i in range(0,Npts): Ans += f(x[i])*w[i] # Sum integrands
print('\n Npts =', Npts, ', Ans =', Ans)
print(' eps =',eps, ', Error =', Ans-(exp(1)-1))
运行结果为
Npts = 10 , Ans = 2.3504023872875797
eps = 3e-14 , Error = 0.6321205588285346
2 逐行解析
Npts = 10; Ans = 0; a = 0.; b = 1.; eps = 3.E-14
w = zeros(2001, float); x = zeros(2001, float) # Arrays
- Npts, Ans, a, b, eps: 初始设置各种变量
- w, x: 一维2001元的浮点数数组，用0.0填充
补充：没想太明白为什么要2001元，感觉元素数量应该和Npts数量一致更好。
def f(x): return exp(x) # Integrand
GaussPoints(Npts, a, b, x, w, eps) # eps: precison of pts
for i in range(0,Npts): Ans += f(x[i])*w[i] # Sum integrands
- f(x): e x e^x ex
- 调用
GaussPoints函数，里面将会把结果保存到x, w这两个数组里面。 - 求和，对应公式为 A n s = ∑ i = 0 N f ( x i ) w i Ans = \sum_{i=0}^{N} f(x_i)w_i Ans=∑i=0Nf(xi)wi
print('\n Npts =', Npts, ', Ans =', Ans)
print(' eps =',eps, ', Error =', Ans-(exp(1)-1))
对计算结果进行输出
其中Error计算对应的公式为
A
n
s
−
(
e
−
1
)
Ans-(e-1)
Ans−(e−1)
1180
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
