python一对一辅导教程：Computational Problems for Physics chapter 1-C Code Listings 1.12 - 1.17

作者自我介绍：大爽歌, 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=hf(t,y_n)$$
$$k_2=hf(t+\frac{h}{2},y_n+\frac{k_1}{2})$$
$$k_3=hf(t+\frac{h}{2},y_n+\frac{k_2}{2})$$
$$k_4=hf(t+\frac{h}{2},y_n+\frac{k_3}{2})$$
$$y_{n+1}=y_n+\frac{k_1}{6}+\frac{k_2}{3}+\frac{k_3}{3}+\frac{k_4}{6}+O(h^5)$$
$$其中O(h^5)代表偏差量，使用公式时不计算，好像也算不了$$

3 逐行解析

解析如下图所示

在每一个阶段，ydumb起的是阶段内部存储暂时值的作用，

4 运行示例

以后面将要调用rk4Algor方法的rk4Call.py中的一次调用的参数为例

• t: 0.0
• h: 0.1
• y: array[ 3. -5.]
• f: $f(t,y)=[y_1,-100y_0-2y_1+10sin\ast (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$
• 第二个元素为 $-100y_0-2y_1+10sin\ast (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)$
• 曲线2添加点 $(t,y_0)$

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,-2g{y}_1-a{y}_0-b{y}_0^3+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)=5si{n}^2(8x)e^{-x^2}-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    

定义函数，内部执行的代码对应公式如下
$$\begin{array}{rl}h& =\frac{B-A}{N-1}\\ strapezoid(A,B,N)& =h(\frac{f(A)+f(B)}{2}+\sum _{i=1}^{N-1}f(A+ih))\end{array}$$

该公式为梯形积分公式，参考文档：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=cos(\pi \frac{i-\frac{1}{4}}{n+\frac{1}{2}})$$

来自直播间老哥的补充：猜测是和勒让得(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\frac{t{p}_N-p_{N-1}}{t^2-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个值，都取值$\frac{2}{(1-t^2)p{p}^2}$和$\frac{2}{(1-t^2)p{p}^2}$

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$
• 调用GaussPoints函数，里面将会把结果保存到x, w这两个数组里面。
• 求和，对应公式为 $Ans={\sum }_{i=0}^{N}f(x_i)w_i$

print('\n Npts =', Npts, ',   Ans =', Ans)
print(' eps =',eps, ', Error =', Ans-(exp(1)-1))

对计算结果进行输出
其中Error计算对应的公式为$Ans-(e-1)$