國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Appendix A
Code use in paper
A.1 Code 01
In[1]:
1 %pylab inline Some Tools
Heun’s Method In[2]:
1 def ode_Heun(f,t0,tf,y0,y1,n):
2 t = np.linspace(t0,tf,n)
3 y = list([y0,y1])
4 for i in range(1,n-1):
5 h = t[i+1]-t[i]
6 fk = f(t[i],y[-1])
7 fk1 = f(t[i+1],y[-1]+h*fk)
8 y.append(y[-1]+h*(fk+fk1)/2.0)
9 return y RK2 Method
‧
2 t = linspace(t0,tf,n)
3 y = list([y0,y1])
4 for i in range(1,n-1):
5 h = t[i+1]-t[i]
6 k1 = h*f(t[i],y[-1])
7 k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0)
8 y.append(y[-1]+k2)
9 return y RK4 Method In[4]:
1 def ode_RK4(f,t0,tf,y0,y1,n):
2 t = linspace(t0,tf,n)
3 y = list([y0,y1])
4 for i in range(1,n-1):
5 h = t[i+1]-t[i]
6 k1 = h*f(t[i],y[-1])
7 k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0)
8 k3 = h*f(t[i]+h/2.0,y[-1]+k2/2.0)
9 k4 = h*f(t[i]+h,y[-1]+k3)
10 y.append(y[-1]+(k1+2*k2+2*k3+k4)/6)
3 for i in range(len(y_real)):
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
4 SSE_Temp = SSE_Temp + (y_real[i]-y_estimate[i])**2
5 SSE_Temp = sqrt(SSE_Temp)
6 return SSE_Temp Data
In[6]:
1 y =
[0.0,0.1,0.3,0.6,0.8,1,2.5,3,4.5,5.7,7.9,9.2,11.5,12,13.2,14.2,15,5.7,16.8,17.5,18.5,19.25,20.1,20.85,21.5,22.15,22.75,23.15,23.5,24,24.4,24.75,25]
2 cdy =
[0.15,0.25,0.25,0.2,0.85,1.0,1.0,1.35,1.7,1.75,1.8,1.4,0.85,1.1,0.9,0.75,0.9,0.9,0.85,0.875,0.8,0.8,0.7,0.65,0.625,0.5,0.375,0.425,0.45,0.375,0.3]
3 n = len(y)
4 n_cdy = len(cdy)
5 t = range(n) In[7]:
1 xlabel(’t:Time (weeks)’)
2 ylabel(’u(t):Height (cm)’)
3 plot(t,y,marker = ’o’)
1. Finite Difference Method 求y′ =αy2+βy 的最佳 α 和 β 利⽤ central difference method
‧
12 z = linalg.solve(A.T*A,A.T*b.T)
13
14 alpha = float(z[0])
15 beta = float(z[1])
16 z
matrix([[-0.00844749],[ 0.21174495]])
Step.2 接下來利⽤上⾯求出的α 和 β，搭配上數值偏微分⽅法來得到⼀組數據逼近
原始資料
也就是說 Given{(yi, ti)}ni=1，我們要⽤ Heun’s Method、RK2 Method 和 RK4 Method
‧
1 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n)
2 SSE_Heun0 = SSE_finitedifference(y_Heun,y)
3 print SSE_Heun0 55.7242174456
RK2 Method In[11]:
1 y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n)
2 SSE_RK2a = SSE_finitedifference(y_RK2,y)
3 print SSE_RK2a 55.66572951
‧
1 y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n)
2 SSE_RK4a = SSE_finitedifference(y_RK4,y)
3 print SSE_RK4a
6 SSE_0 = SSE_finitedifference(y,y_1)
7 print SSE_0 In[14]:
1 count = 0
2 SSE = SSE_0
3 while count <1000:
4 y_temp = [0.0, 0.0]
5
6 # 對 起 始 點 增 加 擾 動
7 e1 = random.uniform(-0.1,0.1)
8 e2 = random.uniform(-0.1,0.1)
‧
9 y_temp[0] = y_1[0]+e1
10 y_temp[1] = y_1[0]+e2
11
12 # ⽣ 成 新 的 Data
13 for i in range(n-2):
14 y_temp.append(y_temp[-2]+2*alpha*y_temp[-1]**2+2*beta*y_temp [-1])
15
16 # 計 算 新 的 最 ⼩ 平 ⽅ 差
17 SSE_temp = SSE_finitedifference(y,y_temp)
18
19 # ⽐ 較 最 ⼩ 平 ⽅ 差
20 if SSE_temp < SSE:
21 SSE = SSE_temp
22 y_1 = y_temp
28 print y_1[0], y_1[1]
5.659286856
-0.0084474869693 0.211744951986 1.16903645355 1.32198681127
Heun’s Method In[15]:
1 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y_1[0],y_1[1],n)
‧
2 SSE_Heun0 = SSE_finitedifference(y_Heun,y)
3 print SSE_Heun0 30.2074082818
RK2 Method In[16]:
1 y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y_1[0],y_1[1],n)
2 SSE_RK2a = SSE_finitedifference(y_RK2,y)
3 print SSE_RK2a 30.1259883441
RK4 Method In[17]:
1 y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y_1[0],y_1[1],n)
2 SSE_RK4a = SSE_finitedifference(y_RK4,y)
3 print SSE_RK4a
2 while flag < 100:
3 count = 0
4 SSE = SSE_0
5 while count <1000:
6 y_temp = [y_1[0], y_1[1]]
‧
9 e1 = random.uniform(-0.001,0.001)
10 e2 = random.uniform(-0.001,0.001)
11 alpha_temp = alpha + e1
12 beta_temp = beta + e2
13
14 # ⽣ 成 新 的 Data
15 for i in range(n-2):
16 y_temp.append(y_temp[-2]+2*alpha_temp*y_temp[-1]**2+2*
beta_temp*y_temp[-1])
17 # 計 算 新 的 最 ⼩ 平 ⽅ 差
18 SSE_temp = SSE_finitedifference(y,y_temp)
19
20 # ⽐ 較 最 ⼩ 平 ⽅ 差
21 if SSE_temp < SSE:
22 SSE = SSE_temp
23 y_1 = y_temp
24 alpha = alpha_temp
25 beta = beta_temp
26
27 count += 1
28
29 count = 0
30 while count <1000:
31 y_temp = [0.0, 0.0]
32
33 # 對 起 始 點 增 加 擾 動
‧
34 e1 = random.uniform(-0.1,0.1)
35 e2 = random.uniform(-0.1,0.1)
36 y_temp[0] = y_1[0]+e1
37 y_temp[1] = y_1[0]+e2
38
39 # ⽣ 成 新 的 Data
40 for i in range(n-2):
41 y_temp.append(y_temp[-2]+2*alpha*y_temp[-1]**2+2*beta*
y_temp[-1])
42
43 # 計 算 新 的 最 ⼩ 平 ⽅ 差
44 SSE_temp = SSE_finitedifference(y,y_temp)
45
46 # ⽐ 較 最 ⼩ 平 ⽅ 差
47 if SSE_temp < SSE:
48 SSE = SSE_temp
49 y_1 = y_temp
55 print y_1[0], y_1[1]
5.39888930477
-0.00959911519495 0.232405855098 0.990787303889 1.07654767318
In[19]:
‧
1 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n)
2 SSE_Heun0 = SSE_finitedifference(y_Heun,y)
3 print SSE_Heun0 5.6053053617
RK2’s Method In[21]:
1 y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n)
2 SSE_RK2a = SSE_finitedifference(y_RK2,y)
3 print SSE_RK2a 5.58681233538
RK4 Method In[22]:
1 y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n)
2 SSE_RK4a = SSE_finitedifference(y_RK4,y)
3 print SSE_RK4a
‧
2 t = np.linspace(t0,tf,n)
3 y = list([y0])
4 for i in range(n-1):
5 h = t[i+1]-t[i]
6 fk = f(t[i],y[-1])
7 fk1 = f(t[i+1],y[-1]+h*fk)
8 y.append(y[-1]+h*(fk+fk1)/2.0)
2 t = linspace(t0,tf,n)
3 y = list([y0])
4 for i in range(n-1):
5 h = t[i+1]-t[i]
6 k1 = h*f(t[i],y[-1])
‧
7 k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0)
8 y.append(y[-1]+k2)
9 return y RK4 Method In[4]:
1 def ode_RK4(f,t0,tf,y0,n):
2 t = linspace(t0,tf,n)
3 y = list([y0])
4 for i in range(n-1):
5 h = t[i+1]-t[i]
6 k1 = h*f(t[i],y[-1])
7 k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0)
8 k3 = h*f(t[i]+h/2.0,y[-1]+k2/2.0)
9 k4 = h*f(t[i]+h,y[-1]+k3)
10 y.append(y[-1]+(k1+2*k2+2*k3+k4)/6)
3 for j in range(len(y_real)):
4 SSE_Temp = SSE_Temp + (y_estimate[j]-y_real[j])**2
5 SSE_Temp = sqrt(SSE_Temp)
6 return SSE_Temp In[6]:
1 def SSE_numerical(y_estimate,y_real,gap):
2 y_T = []
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
3 SSE_Temp = 0.0
4 for i in range(len(y_real)):
5 y_T.append(y_estimate[i*gap])
6 SSE_Temp = SSE_Temp+(y_T[i]-y_real[i])**2
7 SSE_Temp = sqrt(SSE_Temp)
8 return SSE_Temp Data
In[7]:
1 y =
[0.0,0.1,0.3,0.6,0.8,1,2.5,3,4.5,5.7,7.9,9.2,11.5,12,13.2,14.2,15,15.7,16.8,17.5,18.5,19.25,20.1,20.85,21.5,22.15,22.75,23.15,23.5,24,24.4,24.75,25]
2 cdy =
[0.15,0.25,0.25,0.2,0.85,1.0,1.0,1.35,1.7,1.75,1.8,1.4,0.85,1.1,0.9,0.75,0.9,0.9,0.85,0.875,0.8,0.8,0.7,0.65,0.625,0.5,0.375,0.425,0.45,0.375,0.3]
3 n = len(y)
4 n_cdy = len(cdy)
5 t = range(n) In[8]:
1 xlabel(’t:Time (weeks)’)
2 ylabel(’u(t):Height (cm)’)
3 plot(t,y,marker = ’o’)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
2. 利⽤(y, cdy)來找出 y′ =αy2+βy 的最佳 α 和 β Step1. 利⽤最⼩平⽅法求出最佳α 和 β
In[9]:
1 y_5 =
[0.0,0.1,0.3,0.6,0.8,1,2.5,3,4.5,5.7,7.9,9.2,11.5,12,13.2,14.2,15,15.7,16.8,17.5,18.5,19.25,20.1,20.85,21.5,22.15,22.75,23.15,23.5,24,24.4]
2 xlabel(’u(t):Height (cm)’)
3 ylabel(’u(t):velocity (cm/(weeks))’)
4 plot(y_5,cdy,marker = ’o’,color = ’b’,)
In[10]:
1 # 利 ⽤ 最 ⼩ 平 ⽅ 法 求 出 alpha 和 beta
2 A = zeros((n_cdy,2))
3 b = zeros(n_cdy)
‧
5 for i in range(n_cdy):
6 A[i,0] = y[i]**2
13 z = linalg.solve(A.T*A,A.T*b.T)
14
15 alpha = float(z[0])
16 beta = float(z[1])
17 z
5 cdy_estimate = []
6 for i in range(n_cdy):
7 cdy_estimate.append(cdy_fcn(t,y_5[i]))
8
9 xlabel(’u(t):Height (cm)’)
10 ylabel(’u(t):velocity (cm/(weeks))’)
11 plot(y_5,cdy,marker = ’o’,color = ’b’,)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
12 plot(y_5,cdy_estimate,marker = ’o’,color = ’r’,)
In[12]:
1 SSE_cdy0 = SSE_cdy(cdy_estimate,cdy)
2 print SSE_cdy0
3 print alpha, beta 2.31538926492
-0.00863604893773 0.212157561965
Step2. 利⽤求得的 α 和 β，搭配數值微分⽅程的⽅法來找到⼀組函數資料來逼近原
本 Data Setting up In[13]:
1 t0 = 0
2 tf = 30
3 y0 = 1
4
5 # 這 裡 的 gap 必 須 要 是 整 數 ， 代 表 每 格 資 料 間 細 分 的 次 數
6 gap = 100
7 n = 30*gap+1
8 t = linspace(t0,tf,n)
‧
2 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_Heun,n)
3 SSE_Heun0 = SSE_numerical(y_Heun,y_5,gap)
4 print SSE_Heun0 6.3916710696
RK2 Method In[15]:
1 y0_RK2 = y0
2 y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y0_RK2,n)
3 SSE_RK2a = SSE_numerical(y_RK2,y_5,gap)
4 print SSE_RK2a 6.39166483931
RK4 Method In[16]:
1 y0_RK4 = y0
2 y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y0_RK4,n)
3 SSE_RK4a = SSE_numerical(y_RK4,y_5,gap)
4 print SSE_RK4a 6.39165419709
Step.3 對y0做擾動，希望能降低誤差 Heun’s Method
In[17]:
‧
2 SSE_Heun = SSE_Heun0
3 while count <1000:
4
5 # 對 起 始 點 增 加 擾 動
6 e1 = random.uniform(-0.1,0.1)
7 y0_temp = y0_Heun+e1
8
9 # ⽣ 成 新 的 Data
10 y_temp = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_temp,n)
11
12 # 計 算 新 的 最 ⼩ 平 ⽅ 差
13 SSE_temp = SSE_numerical(y_temp,y_5,gap)
14
15 # ⽐ 較 最 ⼩ 平 ⽅ 差
16 if SSE_temp < SSE_Heun:
17 SSE_Heun = SSE_temp
18 y0_Heun = y0_temp
19
‧
2 SSE_RK2 = SSE_RK2a
3 while count <1000:
4
5 # 對 起 始 點 增 加 擾 動
6 e1 = random.uniform(-0.1,0.1)
7 y0_temp = y0_RK2+e1
8
9 # ⽣ 成 新 的 Data
10 y_temp = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_temp,n)
11
12 # 計 算 新 的 最 ⼩ 平 ⽅ 差
13 SSE_temp = SSE_numerical(y_temp,y_5,gap)
14
15 # ⽐ 較 最 ⼩ 平 ⽅ 差
16 if SSE_temp < SSE_RK2:
17 SSE_RK2 = SSE_temp
18 y0_RK2 = y0_temp
19
‧
2 SSE_RK4 = SSE_RK4a
3 while count <1000:
4
5 # 對 起 始 點 增 加 擾 動
6 e1 = random.uniform(-0.1,0.1)
7 y0_temp = y0_RK4+e1
8
9 # ⽣ 成 新 的 Data
10 y_temp = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_temp,n)
11
12 # 計 算 新 的 最 ⼩ 平 ⽅ 差
13 SSE_temp = SSE_numerical(y_temp,y_5,gap)
14
15 # ⽐ 較 最 ⼩ 平 ⽅ 差
16 if SSE_temp < SSE_RK4:
17 SSE_RK4 = SSE_temp
18 y0_RK4 = y0_temp
19
‧
2 t = np.linspace(t0,tf,n)
3 y = list([y0,y1])
4 for i in range(1,n-1):
5 h = t[i+1]-t[i]
6 fk = f(t[i],y[-1])
7 fk1 = f(t[i+1],y[-1]+h*fk)
8 y.append(y[-1]+h*(fk+fk1)/2.0)
2 t = linspace(t0,tf,n)
3 y = list([y0,y1])
4 for i in range(1,n-1):
5 h = t[i+1]-t[i]
6 k1 = h*f(t[i],y[-1])
7 k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0)
8 y.append(y[-1]+k2)
9 return y
‧
2 t = linspace(t0,tf,n)
3 y = list([y0,y1])
4 for i in range(1,n-1):
5 h = t[i+1]-t[i]
6 k1 = h*f(t[i],y[-1])
7 k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0)
8 k3 = h*f(t[i]+h/2.0,y[-1]+k2/2.0)
9 k4 = h*f(t[i]+h,y[-1]+k3)
10 y.append(y[-1]+(k1+2*k2+2*k3+k4)/6)
3 for i in range(len(y_real)):
4 SSE_Temp = SSE_Temp + (y_real[i]-y_estimate[i])**2
5 SSE_Temp = sqrt(SSE_Temp)
6 return SSE_Temp
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
3 n = len(y)
4 n_cdy = len(cdy)
5 t = range(n) In[7]:
1 xlabel(’t:Time (weeks)’)
2 ylabel(’u(t):Height (cm)’)
3 plot(t,y,marker = ’o’)
3. Finite Difference Method 求y′ =αy2+βy+asin(by+c)的最佳α、β、a、b 和 c 利⽤ central difference method
yk+1−yk−1
tk+1−tk−1 =αy2k+βyk+asin(byk+c)
→yk+2=yk+2αy2k+1+2βyk+1+asin(byk+1+c)
Step1. 給定⼀組α、β、a、b 和 c In[8]:
1 alpha = -0.00844749
2 beta = 0.21174495
3 a = 0.0
‧
8 return alpha*y**2+beta*y +a*sin(b*y+c) Heun’s Method
In[10]:
1 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n)
2 SSE_Heun0 = SSE_finitedifference(y_Heun,y)
3 print SSE_Heun0 55.7242229592
RK2 Method In[11]:
1 y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
2 SSE_RK2a = SSE_finitedifference(y_RK2,y)
3 print SSE_RK2a 55.6657350223
RK4 Method In[12]:
1 y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n)
2 SSE_RK4a = SSE_finitedifference(y_RK4,y)
3 print SSE_RK4a 55.3013536756
Step3. 對y0，y1做擾動，希望能降低 SSE
yk+2 =yk+2αy2k+1+2βyk+1+2asin(byk+1+c) In[13]:
1 y_1 = [0.0 , 0.1]
2 for i in range(n-2):
3 y_1.append(y_1[-2]+2*alpha*y_1[-1]**2+2*beta*y_1[-1]+2*a*sin(b*y_1 [-1]+c))
4
5 # 最 ⼩ 平 ⽅ 差
6 SSE_0 = SSE_finitedifference(y,y_1)
7 print SSE_0 64.1919298183
In[14]:
1 count = 0
‧
3 while count <1000:
4 y_temp = [0.0, 0.0]
5
6 # 對 起 始 點 增 加 擾 動
7 e1 = random.uniform(-0.1,0.1)
8 e2 = random.uniform(-0.1,0.1)
9 y_temp[0] = y_1[0]+e1
10 y_temp[1] = y_1[0]+e2
11
12 # ⽣ 成 新 的 Data
13 for i in range(n-2):
14 y_temp.append(y_temp[-2]+2*alpha*y_temp[-1]**2+2*beta*y_temp [-1]+2*a*sin(b*y_temp[-1]+c))
15
16 # 計 算 新 的 最 ⼩ 平 ⽅ 差
17 SSE_temp = SSE_finitedifference(y,y_temp)
18
19 # ⽐ 較 最 ⼩ 平 ⽅ 差
20 if SSE_temp < SSE:
21 SSE = SSE_temp
22 y_1 = y_temp
28 print y_1[0], y_1[1]
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
5.65827878047
-0.00844749 0.21174495 0.0 0.0 0.0 1.20182153439 1.32129224963
In[15]:
1 xlabel(’t:Time (weeks)’)
2 ylabel(’u(t):Height (cm)’)
3 plot(t,y,marker = ’o’,color = ’b’,)
4 plot(t,y_1,marker = ’o’,color = ’r’,)
Heun’s Method In[16]:
1 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y_1[0],y_1[1],n)
2 SSE_Heun0 = SSE_finitedifference(y_Heun,y)
3 print SSE_Heun0 5.79278605472
RK2 Method In[17]:
1 y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y_1[0],y_1[1],n)
2 SSE_RK2a = SSE_finitedifference(y_RK2,y)
‧
1 y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y_1[0],y_1[1],n)
2 SSE_RK4a = SSE_finitedifference(y_RK4,y)
3 print SSE_RK4a
2 while flag < 100:
3 count = 0
4 SSE = SSE_0
5 while count <100:
6 y_temp = [y_1[0], y_1[1]]
7
8 # 對 起 始 點 增 加 擾 動
9 e1 = random.uniform(-0.001,0.001)
10 e2 = random.uniform(-0.001,0.001)
11 e3 = random.uniform(-0.001,0.001)
12 e4 = random.uniform(-0.001,0.001)
13 e5 = random.uniform(-0.001,0.001)
14
15 alpha_temp = alpha + e1
‧
23 for i in range(n-2):
24 y_temp.append(y_temp[-2]+2*alpha_temp*y_temp[-1]**2+2*
beta_temp*y_temp[-1]+a_temp*sin(b_temp*y_temp[-1]+c_temp ))
25 # 計 算 新 的 最 ⼩ 平 ⽅ 差
26 SSE_temp = SSE_finitedifference(y,y_temp)
27
28 # ⽐ 較 最 ⼩ 平 ⽅ 差
29 if SSE_temp < SSE:
30 SSE = SSE_temp
31 y_1 = y_temp
32 alpha = alpha_temp
33 beta = beta_temp
34 a = a_temp
41 while count <100:
‧
45 e1 = random.uniform(-0.1,0.1)
46 e2 = random.uniform(-0.1,0.1)
47 y_temp[0] = y_1[0]+e1
48 y_temp[1] = y_1[1]+e2
49
50 # ⽣ 成 新 的 Data
51 for i in range(n-2):
52 y_temp.append(y_temp[-2]+2*alpha*y_temp[-1]**2+2*beta*
y_temp[-1]+a*sin(b*y_temp[-1]+c))
53
54 # 計 算 新 的 最 ⼩ 平 ⽅ 差
55 SSE_temp = SSE_finitedifference(y,y_temp)
56
57 # ⽐ 較 最 ⼩ 平 ⽅ 差
58 if SSE_temp < SSE:
59 SSE = SSE_temp
60 y_1 = y_temp
66 print y_1[0], y_1[1]
5.40260286876
-0.00951510827755 0.230182300243 -0.0262259316043
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
-0.0147144948651 -0.00726659819372 1.01627076858 1.10278819877
In[20]:
1 xlabel(’t:Time (weeks)’)
2 ylabel(’u(t):Height (cm)’)
3 plot(t,y,marker = ’o’,color = ’b’,)
4 plot(t,y_1,marker = ’o’,color = ’r’,)
In[21]:
1 t0 = t[0]
2 tf = t[-1]
3 y0 = y_1[0]
4 y1 = y_1[1]
5 # 要 解 的 微 分 ⽅ 程
6 def cdy_fcn(t,y):
7 return alpha*y**2+beta*y +a*sin(b*y +c) Heun’s Method
In[22]:
1 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n)
2 SSE_Heun0 = SSE_finitedifference(y_Heun,y)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
3 print SSE_Heun0 5.58823184026
RK2’s Method In[23]:
1 y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n)
2 SSE_RK2a = SSE_finitedifference(y_RK2,y)
3 print SSE_RK2a 5.56985208097
RK4 Method In[24]:
1 y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n)
2 SSE_RK4a = SSE_finitedifference(y_RK4,y)
3 print SSE_RK4a 5.51999849112
A.4 Code 04
In[1]:
1 %pylab inline Some Tools Heun’s Method In[2]:
1 def ode_Heun(f,t0,tf,y0,n):
2 t = np.linspace(t0,tf,n)
‧
7 fk1 = f(t[i+1],y[-1]+h*fk)
8 y.append(y[-1]+h*(fk+fk1)/2.0)
2 t = linspace(t0,tf,n)
3 y = list([y0])
4 for i in range(n-1):
5 h = t[i+1]-t[i]
6 k1 = h*f(t[i],y[-1])
7 k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0)
8 y.append(y[-1]+k2)
9 return y RK4 Method In[4]:
1 def ode_RK4(f,t0,tf,y0,n):
2 t = linspace(t0,tf,n)
3 y = list([y0])
4 for i in range(n-1):
5 h = t[i+1]-t[i]
6 k1 = h*f(t[i],y[-1])
7 k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0)
‧
8 k3 = h*f(t[i]+h/2.0,y[-1]+k2/2.0)
9 k4 = h*f(t[i]+h,y[-1]+k3)
10 y.append(y[-1]+(k1+2*k2+2*k3+k4)/6)
3 for j in range(len(y_real)):
4 SSE_Temp = SSE_Temp + (y_estimate[j]-y_real[j])**2
5 SSE_Temp = sqrt(SSE_Temp)
6 return SSE_Temp In[6]:
1 def SSE_numerical(y_estimate,y_real,gap):
2 y_T = []
3 SSE_Temp = 0.0
4 for i in range(len(y_real)):
5 y_T.append(y_estimate[i*gap])
6 SSE_Temp = SSE_Temp+(y_T[i]-y_real[i])**2
7 SSE_Temp = sqrt(SSE_Temp)
8 return SSE_Temp
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
[0.15,0.25,0.25,0.2,0.85,1.0,1.0,1.35,1.7,1.75,1.8,1.4,0.85,1.1,0.9,0.75,0.9,0.9,0.85,0.875,0.8,0.8,0.7,0.65,0.625,0.5,0.375,0.425,0.45,0.375,0.3]
3 n = len(y)
4 n_cdy = len(cdy)
5 t_0 = range(n) In[8]:
1 xlabel(’t:Time (weeks)’)
2 ylabel(’u(t):Height (cm)’)
3 plot(t_0,y,marker = ’o’)
4. 利⽤(y, cdy)來找出 y′ =αy2+βy+asin(by+c)的最佳α、β、a、b 和 c In[9]:
1 y_5 =
[0.0,0.1,0.3,0.6,0.8,1,2.5,3,4.5,5.7,7.9,9.2,11.5,12,13.2,14.2,15,15.7,16.8,17.5,18.5,19.25,20.1,20.85,21.5,22.15,22.75,23.15,23.5,24,24.4]
2 xlabel(’u(t):Height (cm)’)
3 ylabel(’u(t):velocity (cm/(weeks))’)
4 plot(y_5,cdy,marker = ’o’,color = ’b’,)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Step1. 給定⼀組α、β、a、b 和 c In[10]:
1 alpha = -0.01
2 beta = 0.2
3 a = -0.5
4 b = -0.3
5 c = -0.1 In[11]:
1 def cdy_fcn(t,y):
2 return alpha*y**2+beta*y+a*sin(b*y+c)
3
4 cdy_estimate = []
5 for i in range(n_cdy):
6 cdy_estimate.append(cdy_fcn(t_0,y_5[i]))
7
8 xlabel(’u(t):Height (cm)’)
9 ylabel(’u(t):velocity (cm/(weeks))’)
10 plot(y_5,cdy,marker = ’o’,color = ’b’,)
11 plot(y_5,cdy_estimate,marker = ’o’,color = ’r’,)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
In[12]:
1 SSE_cdy0 = SSE_cdy(cdy_estimate,cdy)
2 print SSE_cdy0
3 print alpha, beta, a, b, c 3.45827066306
-0.01 0.2 -0.5 -0.3 -0.1
Step2. 利⽤求得的α、β、a、b 和 c，搭配數值微分⽅程的⽅法來找到⼀組函數資料
來逼近原本 Data Setting up In[13]:
1 t0 = 0
2 tf = 30
3 y0 = 1
4
5 # 這 裡 的 gap 必 須 要 是 整 數 ， 代 表 每 格 資 料 間 細 分 的 次 數
6 gap = 100
7 n = 30*gap+1
8 t = linspace(t0,tf,n) Heun’s method
‧
2 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_Heun,n)
3 SSE_Heun0 = SSE_numerical(y_Heun,y_5,gap)
4 print SSE_Heun0 22.7715170995
RK2 Method In[15]:
1 y0_RK2 = y0
2 y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y0_RK2,n)
3 SSE_RK2a = SSE_numerical(y_RK2,y_5,gap)
4 print SSE_RK2a 22.7715218758
RK4 Method In[16]:
1 y0_RK4 = y0
2 y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y0_RK4,n)
3 SSE_RK4a = SSE_numerical(y_RK4,y_5,gap)
4 print SSE_RK4a 22.7715246995
Step.3 對α、β、a、b 和 c 做擾動，降低⼆階函數估計值跟實際值
In[17]:
1 count = 0
2 SSE_0 = SSE_cdy0
‧
3 while count <10000:
4
5 # 對 起 始 點 增 加 擾 動
6 e1 = random.uniform(-0.001,0.001)
7 e2 = random.uniform(-0.001,0.001)
8 e3 = random.uniform(-0.001,0.001)
9 e4 = random.uniform(-0.001,0.001)
10 e5 = random.uniform(-0.001,0.001)
11 alpha_temp = alpha + e1
12 beta_temp = beta + e2
22 cdy_estimate_temp = []
23 for i in range(n_cdy):
24 cdy_estimate_temp.append(cdy_fcn_temp(t,y_5[i]))
25
26 # 計 算 新 的 最 ⼩ 平 ⽅ 差
27 SSE_temp = SSE_cdy(cdy_estimate_temp,cdy)
28
29 # ⽐ 較 最 ⼩ 平 ⽅ 差
30 if SSE_temp < SSE_0:
‧
31 SSE_0 = SSE_temp
32 alpha = alpha_temp
33 beta = beta_temp
34 a = a_temp
-0.0114876114475 0.262698012377 -0.53344134757 -0.318721070024 -0.120528934405
In[18]:
1 def cdy_fcn(t,y):
2 return alpha*y**2+beta*y+a*sin(b*y+c) In[19]:
1 cdy_estimate = []
2 for i in range(n_cdy):
3 cdy_estimate.append(cdy_fcn(t,y_5[i]))
4
5 xlabel(’u(t):Height (cm)’)
6 ylabel(’u(t):velocity (cm/(weeks))’)
7 plot(y_5,cdy,marker = ’o’,color = ’b’,)
8 plot(y_5,cdy_estimate,marker = ’o’,color = ’r’,)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
In[20]:
1 t0 = 0
2 tf = 30
3 y0 = 1
4
5 # 這 裡 的 gap 必 須 要 是 整 數 ， 代 表 每 格 資 料 間 細 分 的 次 數
6 gap = 100
7 n = 30*gap+1
8 t = linspace(t0,tf,n) In[21]:
1 y0_Heun = y0
2 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_Heun,n)
3 SSE_Heun0 = SSE_numerical(y_Heun,y_5,gap)
4 print SSE_Heun0 20.2251287211
In[22]:
1 xlabel(’t:Time (weeks)’)
2 ylabel(’u(t):Height (cm)’)
3 plot(t_0,y,marker = ’o’,color = ’blue’)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
4 plot(t,y_Heun,color = ’red’)
Step.4 對y0做擾動，希望能降低誤差 Heun’s Method
In[23]:
1 count = 0
2 SSE_Heun = SSE_Heun0
3 while count <1000:
4
5 # 對 起 始 點 增 加 擾 動
6 e1 = random.uniform(-0.1,0.1)
7 y0_temp = y0_Heun+e1
8
9 # ⽣ 成 新 的 Data
10 y_temp = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_temp,n)
11
12 # 計 算 新 的 最 ⼩ 平 ⽅ 差
13 SSE_temp = SSE_numerical(y_temp,y_5,gap)
14
15 # ⽐ 較 最 ⼩ 平 ⽅ 差
16 if SSE_temp < SSE_Heun:
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
17 SSE_Heun = SSE_temp
18 y0_Heun = y0_temp
19
20 count += 1
21
22 print SSE_Heun
23 print y0_Heun 1.68479162933 0.0489081076939
In[24]:
1 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_Heun,n)
2 xlabel(’t:Time (weeks)’)
3 ylabel(’u(t):Height (cm)’)
4 plot(t_0,y,marker = ’o’,color = ’blue’)
5 plot(t,y_Heun,color = ’red’)
RK2 Method In[25]:
1 count = 0
2 SSE_RK2 = SSE_RK2a
‧
3 while count <1000:
4
5 # 對 起 始 點 增 加 擾 動
6 e1 = random.uniform(-0.1,0.1)
7 y0_temp = y0_RK2+e1
8
9 # ⽣ 成 新 的 Data
10 y_temp = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_temp,n)
11
12 # 計 算 新 的 最 ⼩ 平 ⽅ 差
13 SSE_temp = SSE_numerical(y_temp,y_5,gap)
14
15 # ⽐ 較 最 ⼩ 平 ⽅ 差
16 if SSE_temp < SSE_RK2:
17 SSE_RK2 = SSE_temp
18 y0_RK2 = y0_temp
19
1 y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y0_RK2,n)
2 xlabel(’t:Time (weeks)’)
3 ylabel(’u(t):Height (cm)’)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
4 plot(t_0,y,marker = ’o’,color = ’blue’)
5 plot(t,y_RK2,color = ’red’)
RK4 Method In[27]:
1 count = 0
2 SSE_RK4 = SSE_RK4a
3 while count <1000:
4
5 # 對 起 始 點 增 加 擾 動
6 e1 = random.uniform(-0.1,0.1)
7 y0_temp = y0_RK4+e1
8
9 # ⽣ 成 新 的 Data
10 y_temp = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_temp,n)
11
12 # 計 算 新 的 最 ⼩ 平 ⽅ 差
13 SSE_temp = SSE_numerical(y_temp,y_5,gap)
14
15 # ⽐ 較 最 ⼩ 平 ⽅ 差
16 if SSE_temp < SSE_RK4:
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
17 SSE_RK4 = SSE_temp
18 y0_RK4 = y0_temp
19
20 count += 1
21
22 print SSE_RK4
23 print y0_RK4 1.68491878659 0.0486656124813
In[28]:
1 y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y0_RK4,n)
2 xlabel(’t:Time (weeks)’)
3 ylabel(’u(t):Height (cm)’)
4 plot(t_0,y,marker = ’o’,color = ’blue’)
5 plot(t,y_RK4,color = ’red’)
A.5 Code 05
In[1]:
1 %pylab inline
‧
2 t = np.linspace(t0,tf,n)
3 y = list([y0])
4 for i in range(n-1):
5 h = t[i+1]-t[i]
6 fk = f(t[i],y[-1])
7 fk1 = f(t[i+1],y[-1]+h*fk)
8 y.append(y[-1]+h*(fk+fk1)/2.0)
2 t = linspace(t0,tf,n)
3 y = list([y0])
4 for i in range(n-1):
5 h = t[i+1]-t[i]
6 k1 = h*f(t[i],y[-1])
7 k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0)
8 y.append(y[-1]+k2)
9 return y RK4 Method In[4]:
1 def ode_RK4(f,t0,tf,y0,n):
2 t = linspace(t0,tf,n)
‧
7 k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0)
8 k3 = h*f(t[i]+h/2.0,y[-1]+k2/2.0)
9 k4 = h*f(t[i]+h,y[-1]+k3)
10 y.append(y[-1]+(k1+2*k2+2*k3+k4)/6)
3 for j in range(len(y_real)):
4 SSE_Temp = SSE_Temp + (y_estimate[j]-y_real[j])**2
5 SSE_Temp = sqrt(SSE_Temp)
6 return SSE_Temp In[6]:
1 def SSE_numerical(y_estimate,y_real,gap):
2 y_T = []
3 SSE_Temp = 0.0
4 for i in range(len(y_real)):
5 y_T.append(y_estimate[i*gap])
6 SSE_Temp = SSE_Temp+(y_T[i]-y_real[i])**2
7 SSE_Temp = sqrt(SSE_Temp)
8 return SSE_Temp Data
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
In[7]:
1 y =
[0.0,0.1,0.3,0.6,0.8,1,2.5,3,4.5,5.7,7.9,9.2,11.5,12,13.2,14.2,15,15.7,16.8,17.5,18.5,19.25,20.1,20.85,21.5,22.15,22.75,23.15,23.5,24,24.4,24.75,25]
2 cdy =
[0.15,0.25,0.25,0.2,0.85,1.0,1.0,1.35,1.7,1.75,1.8,1.4,0.85,1.1,0.9,0.75,0.9,0.9,0.85,0.875,0.8,0.8,0.7,0.65,0.625,0.5,0.375,0.425,0.45,0.375,0.3]
3 n = len(y)
4 n_cdy = len(cdy)
5 t_0 = range(n) In[8]:
1 xlabel(’t:Time (weeks)’)
2 ylabel(’u(t):Height (cm)’)
3 plot(t_0,y,marker = ’o’)
5. 利⽤(y, cdy)來找出 y′ =αy2+βy+asin(by+c)的最佳α、β、a、b 和 c In[9]:
1 y_5 =
[0.0,0.1,0.3,0.6,0.8,1,2.5,3,4.5,5.7,7.9,9.2,11.5,12,13.2,14.2,15,15.7,16.8,17.5,18.5,19.25,20.1,20.85,21.5,22.15,22.75,23.15,23.5,24,24.4]
2 xlabel(’u(t):Height (cm)’)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
3 ylabel(’u(t):velocity (cm/(weeks))’)
4 plot(y_5,cdy,marker = ’o’,color = ’b’,)
Step1. 利⽤ Matlab Curve Fitting Tool 找到⼀組α、β、a、b 和 c In[10]:
1 alpha = -0.011
2 beta = 0.2612
3 a = -0.5341
4 b = -0.3127
5 c = -0.1563 In[11]:
1 def cdy_fcn(t,y):
2 return alpha*y**2+beta*y+a*sin(b*y+c)
3
4 cdy_estimate = []
5 for i in range(n_cdy):
6 cdy_estimate.append(cdy_fcn(t_0,y_5[i]))
7
8 xlabel(’u(t):Height (cm)’)
9 ylabel(’u(t):velocity (cm/(weeks))’)
10 plot(y_5,cdy,marker = ’o’,color = ’b’,)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
11 plot(y_5,cdy_estimate,marker = ’o’,color = ’r’,)
In[12]:
1 SSE_cdy0 = SSE_cdy(cdy_estimate,cdy)
2 print SSE_cdy0
3 print alpha, beta, a, b, c 1.11321312693
-0.011 0.2612 -0.5341 -0.3127 -0.1563
Step2. 利⽤上述的α、β、a、b 和 c，搭配數值微分⽅程的⽅法來找到⼀組函數資料
來逼近原本 Data Setting up In[13]:
1 t0 = 0
2 tf = 30
3 y0 = 1
4
5 # 這 裡 的 gap 必 須 要 是 整 數 ， 代 表 每 格 資 料 間 細 分 的 次 數
6 gap = 100
7 n = 30*gap+1
8 t = linspace(t0,tf,n)
‧
2 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_Heun,n)
3 SSE_Heun0 = SSE_numerical(y_Heun,y_5,gap)
4 print SSE_Heun0 22.60031044
RK2 Method In[15]:
1 y0_RK2 = y0
2 y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y0_RK2,n)
3 SSE_RK2a = SSE_numerical(y_RK2,y_5,gap)
4 print SSE_RK2a 22.6003363704
RK4 Method In[16]:
1 y0_RK4 = y0
2 y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y0_RK4,n)
3 SSE_RK4a = SSE_numerical(y_RK4,y_5,gap)
4 print SSE_RK4a 22.6003603238
Step.4 對y0做擾動，希望能降低誤差 Heun’s Method
In[17]:
‧
2 SSE_Heun = SSE_Heun0
3 while count <1000:
4
5 # 對 起 始 點 增 加 擾 動
6 e1 = random.uniform(-0.1,0.1)
7 y0_temp = y0_Heun+e1
8
9 # ⽣ 成 新 的 Data
10 y_temp = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_temp,n)
11
12 # 計 算 新 的 最 ⼩ 平 ⽅ 差
13 SSE_temp = SSE_numerical(y_temp,y_5,gap)
14
15 # ⽐ 較 最 ⼩ 平 ⽅ 差
16 if SSE_temp < SSE_Heun:
17 SSE_Heun = SSE_temp
18 y0_Heun = y0_temp
19
1 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_Heun,n)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
2 xlabel(’t:Time (weeks)’)
3 ylabel(’u(t):Height (cm)’)
4 plot(t_0,y,marker = ’o’,color = ’blue’)
5 plot(t,y_Heun,color = ’red’)
RK2 Method In[19]:
1 count = 0
2 SSE_RK2 = SSE_RK2a
3 while count <1000:
4
5 # 對 起 始 點 增 加 擾 動
6 e1 = random.uniform(-0.1,0.1)
7 y0_temp = y0_RK2+e1
8
9 # ⽣ 成 新 的 Data
10 y_temp = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_temp,n)
11
12 # 計 算 新 的 最 ⼩ 平 ⽅ 差
13 SSE_temp = SSE_numerical(y_temp,y_5,gap)
14
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
15 # ⽐ 較 最 ⼩ 平 ⽅ 差
16 if SSE_temp < SSE_RK2:
17 SSE_RK2 = SSE_temp
18 y0_RK2 = y0_temp
19
20 count += 1
21
22 print SSE_RK2
23 print y0_RK2 2.12599915696 -0.0153907948845
In[20]:
1 y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y0_RK2,n)
2 xlabel(’t:Time (weeks)’)
3 ylabel(’u(t):Height (cm)’)
4 plot(t_0,y,marker = ’o’,color = ’blue’)
5 plot(t,y_RK2,color = ’red’)
RK4 Method In[21]:
‧
2 SSE_RK4 = SSE_RK4a
3 while count <1000:
4
5 # 對 起 始 點 增 加 擾 動
6 e1 = random.uniform(-0.1,0.1)
7 y0_temp = y0_RK4+e1
8
9 # ⽣ 成 新 的 Data
10 y_temp = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0_temp,n)
11
12 # 計 算 新 的 最 ⼩ 平 ⽅ 差
13 SSE_temp = SSE_numerical(y_temp,y_5,gap)
14
15 # ⽐ 較 最 ⼩ 平 ⽅ 差
16 if SSE_temp < SSE_RK4:
17 SSE_RK4 = SSE_temp
18 y0_RK4 = y0_temp
19
1 y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y0_RK4,n)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
2 xlabel(’t:Time (weeks)’)
3 ylabel(’u(t):Height (cm)’)
4 plot(t_0,y,marker = ’o’,color = ’blue’)
5 plot(t,y_RK4,color = ’red’)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
Bibliography
[1] 海岸綠提—⽔筆仔 http://163.20.52.80/stu635/cwpspage/mang/study/index.htm.
[2] William E Boyce. Elementary differential equations and boundary value problems. Wiley, 9th edition, 2010.
[3] Brian Bradie. A friendly introduction to nummerical analysis. Pearson Education, Inc., 2006.
[4] Laurence D Hoffmann. Applied calculus for bussiness, economics, and the social and life sciences. McGraw-Hill, eleventh edition, 2014.
[5] Tzeng Jeng-nan. 數 值 微 分 http://glophy.com/index.php/2014-02-07-01-06-58/2014-02-07-01-07-46/79-2014-02-05-07-28-44.
[6] David Kincaid. Numerical analysis: Mathematics of scientific computing. Brooks/Cole, third edition, 2002.
[7] Chen Ren-fa. Nonlinear differential equation of second order and its applications. 2015.