数値積分
モジュールのインポート
import numpy as np
from scipy import integrate # 数値積分
import matplotlib.pyplot as plt
import scienceplots
#Warning : As of version 2.0.0, you need to add import scienceplots before setting the style (plt.style.use('science')).
plt.style.use(['science', 'notebook'])
# rcPramsを使ってTeXのフォントをデフォルトにする
#plt.rcParams['font.family'] = 'Times New Roman' # font familyの設定
plt.rcParams['font.family'] = 'serif' # font familyの設定
plt.rcParams['mathtext.fontset'] = 'cm' # math fontの設定
混合型台形則
def daikei(ff, n, a, b): # 混合型台形則
h = (b-a)/(n-1)
y = np.sum(ff[0:n])
y = y-(ff[0]+ff[n-1])*0.5
y = h*y
return y
混合型シンプソン1/3則
def simpson(ff, n, a, b): # 混合型シンプソン1/3則
h = (b-a)/(n-1)
if n==2:# 台形則
y = (ff[0]+ff[n-1])*h*0.5
elif np.mod(n,2)==1:# 混合型シンプソンの1/3則
y4 = np.sum(ff[1:n:2])
y2 = np.sum(ff[2:n-1:2])
y = ff[0]+4.0*y4+2.0*y2+ff[n-1]
y = y*h/3.0
else:# 混合型シンプソンの1/3則 + 3/8則
y = ff[0]+3.0*(ff[1]+ff[2])+ff[3]
y = y*h*3.0/8.0
if n>=6:
y4 = np.sum(ff[4:n:2])
y2 = np.sum(ff[5:n-1:2])
y = y +(ff[3]+4.0*y4+2.0*y2+ff[n-1])*h/3.0
return y
被積分関数
関数 $y = \sqrt{4-x^2}$ を定義して，
def func3(x): # 1変数の被積分関数
y = np.sqrt(4.0-x**2)
return y
積分範囲を$[0,2]$として，
q = 180.0/np.pi
a = 0.0
b = 2.0
nx = 101
print("f(x) = sqrt(4-x**2)")
print(f"{'n':>4s} {'daikei':>12s} {'error1':>12s} {'simpson':>12s} {'error2':>12s}")
for i in range(2,nx):
n = i
h = (b-a)/(n-1)
xx = np.linspace(a,b,n)
ff = np.empty(n)
ff[:] = [func3(xx[i]) for i in range(n)]
y1 = daikei(ff, n, a, b)
y2 = simpson(ff, n, a, b)
print(f"{n:>4.0f} {y1:>12.9f} {y1-np.pi:>12.9f} {y2:>12.9f} {y2-np.pi:>12.9f}")
f(x) = sqrt(4-x**2)
n daikei error1 simpson error2
2 2.000000000 -1.141592654 2.000000000 -1.141592654
3 2.732050808 -0.409541846 2.976067743 -0.165524910
4 2.917553379 -0.224039275 3.032247551 -0.109345102
5 2.995709068 -0.145883585 3.083595155 -0.057997499
6 3.037048829 -0.104543825 3.100013266 -0.041579388
7 3.061983022 -0.079609631 3.110126237 -0.031466417
8 3.078371909 -0.063220744 3.116625004 -0.024967650
9 3.089819144 -0.051773509 3.121189170 -0.020403484
10 3.098185383 -0.043407271 3.124497987 -0.017094667
11 3.104518326 -0.037074327 3.127008159 -0.014584495
12 3.109448305 -0.032144348 3.128954030 -0.012638624
13 3.113374919 -0.028217734 3.130505552 -0.011087102
14 3.116562476 -0.025030177 3.131761799 -0.009830855
15 3.119192049 -0.022400605 3.132798762 -0.008793892
16 3.121391412 -0.020201241 3.133664561 -0.007928093
17 3.123253038 -0.018339616 3.134397669 -0.007194985
18 3.124845323 -0.016747331 3.135023924 -0.006568730
19 3.126219832 -0.015372822 3.135564648 -0.006028005
20 3.127416097 -0.014176556 3.136034808 -0.005557845
21 3.128464880 -0.013127774 3.136447064 -0.005145589
22 3.129390437 -0.012202217 3.136810615 -0.004782038
23 3.130212126 -0.011380528 3.137133399 -0.004459254
24 3.130945564 -0.010647090 3.137421346 -0.004171307
25 3.131603473 -0.009989181 3.137679658 -0.003912996
26 3.132196311 -0.009396343 3.137912309 -0.003680345
27 3.132732744 -0.008859909 3.138122834 -0.003469820
28 3.133220011 -0.008372643 3.138313987 -0.003278667
29 3.133664195 -0.007928458 3.138488244 -0.003104409
30 3.134070449 -0.007522205 3.138647567 -0.002945086
31 3.134443156 -0.007149498 3.138793737 -0.002798917
32 3.134786070 -0.006806583 3.138928184 -0.002664470
33 3.135102423 -0.006490231 3.139052218 -0.002540436
34 3.135395007 -0.006197646 3.139166905 -0.002425749
35 3.135666251 -0.005926403 3.139273227 -0.002319427
36 3.135918271 -0.005674383 3.139371993 -0.002220661
37 3.136152922 -0.005439732 3.139463952 -0.002128702
38 3.136371835 -0.005220819 3.139549727 -0.002042927
39 3.136576449 -0.005016205 3.139629899 -0.001962755
40 3.136768038 -0.004824616 3.139704955 -0.001887698
41 3.136947734 -0.004644919 3.139775352 -0.001817301
42 3.137116546 -0.004476107 3.139841476 -0.001751178
43 3.137275375 -0.004317278 3.139903688 -0.001688966
44 3.137425028 -0.004167625 3.139962299 -0.001630355
45 3.137566231 -0.004026423 3.140017599 -0.001575054
46 3.137699636 -0.003893017 3.140069840 -0.001522814
47 3.137825834 -0.003766819 3.140119258 -0.001473396
48 3.137945359 -0.003647295 3.140166057 -0.001426596
49 3.138058693 -0.003533961 3.140210433 -0.001382221
50 3.138166277 -0.003426377 3.140252553 -0.001340101
51 3.138268511 -0.003324142 3.140292578 -0.001300076
52 3.138365760 -0.003226893 3.140330648 -0.001262006
53 3.138458359 -0.003134295 3.140366897 -0.001225757
54 3.138546612 -0.003046042 3.140401442 -0.001191212
55 3.138630799 -0.002961854 3.140434395 -0.001158258
56 3.138711179 -0.002881474 3.140465856 -0.001126797
57 3.138787988 -0.002804666 3.140495919 -0.001096735
58 3.138861444 -0.002731209 3.140524667 -0.001067986
59 3.138931749 -0.002660905 3.140552182 -0.001040471
60 3.138999088 -0.002593565 3.140578535 -0.001014119
61 3.139063635 -0.002529019 3.140603794 -0.000988859
62 3.139125547 -0.002467106 3.140628022 -0.000964632
63 3.139184975 -0.002407679 3.140651276 -0.000941377
64 3.139242054 -0.002350600 3.140673611 -0.000919043
65 3.139296913 -0.002295741 3.140695076 -0.000897578
66 3.139349671 -0.002242983 3.140715718 -0.000876935
67 3.139400438 -0.002192215 3.140735582 -0.000857072
68 3.139449319 -0.002143334 3.140754706 -0.000837947
69 3.139496410 -0.002096243 3.140773130 -0.000819523
70 3.139541802 -0.002050852 3.140790888 -0.000801765
71 3.139585579 -0.002007075 3.140808015 -0.000784639
72 3.139627820 -0.001964834 3.140824540 -0.000768114
73 3.139668600 -0.001924054 3.140840493 -0.000752161
74 3.139707989 -0.001884665 3.140855901 -0.000736753
75 3.139746052 -0.001846602 3.140870790 -0.000721863
76 3.139782850 -0.001809803 3.140885185 -0.000707469
77 3.139818443 -0.001774211 3.140899107 -0.000693546
78 3.139852884 -0.001739770 3.140912579 -0.000680075
79 3.139886225 -0.001706429 3.140925620 -0.000667033
80 3.139918514 -0.001674139 3.140938250 -0.000654404
81 3.139949798 -0.001642856 3.140950486 -0.000642168
82 3.139980120 -0.001612534 3.140962345 -0.000630308
83 3.140009520 -0.001583134 3.140973844 -0.000618809
84 3.140038037 -0.001554616 3.140984997 -0.000607656
85 3.140065709 -0.001526945 3.140995820 -0.000596834
86 3.140092569 -0.001500085 3.141006325 -0.000586329
87 3.140118651 -0.001474003 3.141016525 -0.000576129
88 3.140143986 -0.001448668 3.141026433 -0.000566221
89 3.140168603 -0.001424051 3.141036060 -0.000556594
90 3.140192531 -0.001400123 3.141045417 -0.000547236
91 3.140215796 -0.001376858 3.141054516 -0.000538138
92 3.140238424 -0.001354230 3.141063364 -0.000529289
93 3.140260439 -0.001332215 3.141071973 -0.000520680
94 3.140281863 -0.001310790 3.141080351 -0.000512302
95 3.140302720 -0.001289934 3.141088507 -0.000504147
96 3.140323029 -0.001269624 3.141096449 -0.000496205
97 3.140342811 -0.001249842 3.141104184 -0.000488470
98 3.140362085 -0.001230569 3.141111720 -0.000480933
99 3.140380868 -0.001211786 3.141119065 -0.000473589
100 3.140399178 -0.001193475 3.141126224 -0.000466429
Numpy と Scipy の台形公式の数値積分を比較すると，
y4n = np.trapz(ff, xx) # Numpyの台形公式 (trapezoidal rule)
print("np.trapz: ",y4n,y4n-np.pi) # 積分値と誤差
np.trapz: 3.140399178114616 -0.0011934754751772303
https://docs.scipy.org/doc//scipy-1.4.1/reference/generated/scipy.integrate.trapz.html
y4 = integrate.trapz(ff, xx) # Scipyの台形公式 (trapezoidal rule)
print("integrate.trapz: ",y4,y4-np.pi) # 積分値と誤差
integrate.trapz: 3.140399178114616 -0.0011934754751772303
Scipy のシンプソンの公式の数値積分より，
y5 = integrate.simps(ff, xx) # シンプソンの公式 (Simpson's rule)
print("integrate.simps: ",y5,y5-np.pi) # 積分値と誤差
integrate.simps: 3.140876361334483 -0.0007162922553103357
Fortran のライブラリ QUADPACK を使ったガウスの数値積分 (Gaussian quadrature)より，
y3 = integrate.quad(func3, a, b) # 1次元積分(適応積分)
print("integrate.quad :",y3) # 積分値と推定誤差
integrate.quad : (3.1415926535897922, 3.533564552071766e-10)
fig = plt.figure() # グラフ領域の作成
plt.plot(xx,ff)
plt.axis('square')
fig.savefig('p2_ex03_1.pdf')
2重積分
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.dblquad.html
def func33(y,x): # 2変数の被積分関数
f = np.sqrt(x**2+y**2)
return f
x_min, x_max = 0.0, 1.0
y33 = integrate.dblquad(func33, x_min, x_max, 0.0, lambda x: np.sqrt(1-x**2)) # 2次元積分
print("integrate.dblquad:",y33) # 積分値と推定誤差
print("sol:",1.0/6.0*np.pi)
integrate.dblquad: (0.5235987755981053, 1.4792290955080084e-08)
sol: 0.5235987755982988
n = 101
x = np.linspace(-0.5, 1.5, n)
y = np.linspace(-1, 1, n)
X,Y = np.meshgrid(x, y)
v = np.linspace(0, 2, 201)
fig = plt.figure()
plt.xlabel(r"$x$") # x軸のラベル設定
plt.ylabel(r"$y$") # y軸のラベル設定
Cf = plt.contourf(X, Y, func33(Y, X), 8, levels=v, alpha=1.0, cmap="rainbow") # plt.cm.hot
C = plt.contour(X, Y, func33(Y, X), 8, colors='gray', linewidths=0.8)
cbar = plt.colorbar(Cf, pad=0.1, shrink=1.0, orientation="vertical", ticks=[0,0.5, 1, 1.5, 2])
#cbar = plt.colorbar(Cf, pad=0.15, shrink=0.6, orientation="horizontal", ticks=[0,0.5, 1, 1.5, 2])
plt.clabel(C, inline=1, fontsize=12)
plt.axis('square')
fig.savefig('p2_ex03_2.pdf')
plt.show()
前のページに戻る
「目次」のページに戻る