使用Python来解决一个欧拉法画图的问题

Euler's Method with Python

Exercise 1.1 (Page 15)

The velocity of a freely falling ball near Earth's surface is described by the equation

$$\frac{\mathrm{d}v}{\mathrm{d}t} = -g$$

calculate $v$ as a function of $t$

For simplicity, assume that the initial velocity is zero, and calculate the solution for times $t=0$ to $t=10$ s

Euler's Method in Problem 1.1

方程和初始条件为

$$\frac{\mathrm{d}v}{\mathrm{d}t} = -g, v(0) = 0$$

对$t_0 = 0$附近的$t_1$，有：

$$v(t_1) = v'(t_0)(t_1 - t_0) + v_0 = -g \Delta t + v_0$$

同理，对于后续的$t_j$，均有

$$v(t_{j+1}) = -g \Delta t + v_j$$

由此计算出$t_1 , v_1$、$t_2 , v_2$……$t_j , v_j$等一系列点组成的曲线可近似视作$v$随$t$变化的函数图像。

Solution with Python

import numpy as np
from matplotlib import pyplot as plt

引入完了以后给它一个初始值，取$t = 0$到$t = 10$之间的100个点：

t0 = 0
v0 = 0
tf = 10.0
n = 101
g = 9.8
deltat = (tf - t0)/(n-1)

设置坐标轴：

t = np.linspace(t0,tf,n)

v = np.zeros([n])

用欧拉法循环算点：

v[0] = v0

for j in range(1,n):
    v[j] = (-g*deltat) + v[j-1]

for j in range(1,n):
    print(t[j],v[j])

(0.10000000000000001, -0.98000000000000009)
(0.20000000000000001, -1.9600000000000002)
(0.30000000000000004, -2.9400000000000004)
(0.40000000000000002, -3.9200000000000004)
(0.5, -4.9000000000000004)
(0.60000000000000009, -5.8800000000000008)
(0.70000000000000007, -6.8600000000000012)
(0.80000000000000004, -7.8400000000000016)
(0.90000000000000002, -8.8200000000000021)
(1.0, -9.8000000000000025)
(1.1000000000000001, -10.780000000000003)
(1.2000000000000002, -11.760000000000003)
(1.3, -12.740000000000004)
(1.4000000000000001, -13.720000000000004)
(1.5, -14.700000000000005)
(1.6000000000000001, -15.680000000000005)
(1.7000000000000002, -16.660000000000004)
(1.8, -17.640000000000004)
(1.9000000000000001, -18.620000000000005)
(2.0, -19.600000000000005)
(2.1000000000000001, -20.580000000000005)
(2.2000000000000002, -21.560000000000006)
(2.3000000000000003, -22.540000000000006)
(2.4000000000000004, -23.520000000000007)
(2.5, -24.500000000000007)
(2.6000000000000001, -25.480000000000008)
(2.7000000000000002, -26.460000000000008)
(2.8000000000000003, -27.440000000000008)
(2.9000000000000004, -28.420000000000009)
(3.0, -29.400000000000009)
(3.1000000000000001, -30.38000000000001)
(3.2000000000000002, -31.36000000000001)
(3.3000000000000003, -32.340000000000011)
(3.4000000000000004, -33.320000000000007)
(3.5, -34.300000000000004)
(3.6000000000000001, -35.280000000000001)
(3.7000000000000002, -36.259999999999998)
(3.8000000000000003, -37.239999999999995)
(3.9000000000000004, -38.219999999999992)
(4.0, -39.199999999999989)
(4.1000000000000005, -40.179999999999986)
(4.2000000000000002, -41.159999999999982)
(4.2999999999999998, -42.139999999999979)
(4.4000000000000004, -43.119999999999976)
(4.5, -44.099999999999973)
(4.6000000000000005, -45.07999999999997)
(4.7000000000000002, -46.059999999999967)
(4.8000000000000007, -47.039999999999964)
(4.9000000000000004, -48.01999999999996)
(5.0, -48.999999999999957)
(5.1000000000000005, -49.979999999999954)
(5.2000000000000002, -50.959999999999951)
(5.3000000000000007, -51.939999999999948)
(5.4000000000000004, -52.919999999999945)
(5.5, -53.899999999999942)
(5.6000000000000005, -54.879999999999939)
(5.7000000000000002, -55.859999999999935)
(5.8000000000000007, -56.839999999999932)
(5.9000000000000004, -57.819999999999929)
(6.0, -58.799999999999926)
(6.1000000000000005, -59.779999999999923)
(6.2000000000000002, -60.75999999999992)
(6.3000000000000007, -61.739999999999917)
(6.4000000000000004, -62.719999999999914)
(6.5, -63.69999999999991)
(6.6000000000000005, -64.679999999999907)
(6.7000000000000002, -65.659999999999911)
(6.8000000000000007, -66.639999999999915)
(6.9000000000000004, -67.619999999999919)
(7.0, -68.599999999999923)
(7.1000000000000005, -69.579999999999927)
(7.2000000000000002, -70.559999999999931)
(7.3000000000000007, -71.539999999999935)
(7.4000000000000004, -72.519999999999939)
(7.5, -73.499999999999943)
(7.6000000000000005, -74.479999999999947)
(7.7000000000000002, -75.459999999999951)
(7.8000000000000007, -76.439999999999955)
(7.9000000000000004, -77.419999999999959)
(8.0, -78.399999999999963)
(8.0999999999999996, -79.379999999999967)
(8.2000000000000011, -80.359999999999971)
(8.3000000000000007, -81.339999999999975)
(8.4000000000000004, -82.319999999999979)
(8.5, -83.299999999999983)
(8.5999999999999996, -84.279999999999987)
(8.7000000000000011, -85.259999999999991)
(8.8000000000000007, -86.239999999999995)
(8.9000000000000004, -87.219999999999999)
(9.0, -88.200000000000003)
(9.0999999999999996, -89.180000000000007)
(9.2000000000000011, -90.160000000000011)
(9.3000000000000007, -91.140000000000015)
(9.4000000000000004, -92.120000000000019)
(9.5, -93.100000000000023)
(9.6000000000000014, -94.080000000000027)
(9.7000000000000011, -95.060000000000031)
(9.8000000000000007, -96.040000000000035)
(9.9000000000000004, -97.020000000000039)
(10.0, -98.000000000000043)

输出图像：

plt.plot(t,v,'o')
plt.xlabel("Value of t")
plt.ylabel("Value of v")
plt.title("Approximate Solution with Euler's Method in Problem 1.1")
plt.show()

Reference material and Problems Encountered

总体思路和代码参考的是这位大兄弟的文档，但是在改写运行的过程中遇到了几个小问题。

浮点数，啊Py那坑爹的浮点数

In[48]的初始值赋值，一开始是这么写的：

t0 = 0
v0 = 0
tf = 10
n = 101
g = 9.8
deltat = (tf - t0)/(n-1)

其实看起来蛮对的（符合人类常识的对），然后run的时候就出问题了，我的v[j]出来全是0……
很奇怪啊，我公式应该没写错才对啊，反复看公式：

v[0] = v0

for j in range(1,n):
    v[j] = (-g*deltat) + v[j-1]

就算没有直接跟t[j]有关那v[j]也不至于全是0啊……
然后在跟程序媛小伙伴讨论了半天以后排除了一开始的猜想”那个v[j] 没有取到刷新值，一直取得default“，接着，某人突然福至心灵地想起了程序好像是直接舍去位数的，莫不是$\Delta t$太小了直接被舍成0了？

……$10/100 = 0.1$，哪里太小了，不过这个思路是对的，拿去print了一下，果然deltat是一个int的0，于是把In[48]的tf改写了一下：

t0 = 0
v0 = 0
tf = 10.0
n = 101
g = 9.8
deltat = (tf - t0)/(n-1)

很好，deltat终于变成一个浮点数了，皆大欢喜皆大欢喜。

编码问题

这里问题出在Plot的部分，最初的代码是这样的：

plt.plot(t,v,'o')
plt.xlabel("Value of t")
plt.ylabel("Value of v")
plt.title("Approximate Solution with Euler’s Method in Problem 1.1")
plt.show()

……看起来跟后来的代码没啥差别对不对。当时的报错是：

'ascii' codec can't decode byte 0xe2 in position 31: ordinal not in range(128)

指向上面的第四行，看到ASCII大概知道应该是编码问题，但我这不是英文标题吗怎么也有这个问题的？

哦最后发现那个’s的’打成了中文的标点，改成'就没事儿了。某种意义上考验眼神……

如果要给图的标题或者坐标标签写中文的话，查了一下，按这篇文章的方法在文档开头声明一下非ASCII字符就可以了。