# 試験解答例18¶

## 簡単な行列計算¶

In [1]:
import numpy as np
from pprint import pprint
aa = np.array([[4,-1,-1], [1,2,-1],[3,-1,0]])
pprint(aa)

array([[ 4, -1, -1],
[ 1,  2, -1],
[ 3, -1,  0]])

In [2]:
import scipy.linalg as linalg
inv_a = linalg.inv(aa)
pprint(inv_a)

array([[-0.16666667,  0.16666667,  0.5       ],
[-0.5       ,  0.5       ,  0.5       ],
[-1.16666667,  0.16666667,  1.5       ]])

In [3]:
pprint(np.dot(inv_a,aa))

array([[ 1.,  0.,  0.],
[ 0.,  1.,  0.],
[ 0.,  0.,  1.]])


## 数値解の収束性¶

In [4]:
import numpy as np

def func(x):
return -4*np.exp(-x)+2*np.exp(-2*x)

def df(x):
return 4*np.exp(-x) - 4*np.exp(-2*x)

from scipy.optimize import fsolve
x0 = fsolve(func, -1.0)[0]
print(x0)

-0.69314718056

In [5]:
import matplotlib.pyplot as plt

x1=-1.0
x2=1.0
x = np.linspace(x1, x2, 100)
y = func(x)
plt.plot(x, y, color = 'k')
plt.plot([x1,x2],[0,0])
plt.grid()
plt.show()

<matplotlib.figure.Figure at 0x10a237048>
In [6]:
x1, x2 = -1.0, 0.0
f1, f2 = func(x1), func(x2)
print('%+15s %+15s %+15s %+15s'  % ('x1','x2','f1','f2'))
print('%+15.10f %+15.10f %+15.10f %+15.10f' % (x1,x2,f1,f2))

list_bisec = [[0],[abs(x1-x0)]]
for i in range(0, 10):
x = (x1 + x2)/2
f = func(x)
if (f*f1>=0.0):
x1, f1 = x, f
list_bisec[0].append(i)
list_bisec[1].append(abs(x1-x0))
else:
x2, f2 = x, f
list_bisec[0].append(i)
list_bisec[1].append(abs(x2-x0))

print('%+15.10f %+15.10f %+15.10f %+15.10f' % (x1,x2,f1,f2))

print(list_bisec)

             x1              x2              f1              f2
-1.0000000000   +0.0000000000   +3.9049848840   -2.0000000000
-1.0000000000   -0.5000000000   +3.9049848840   -1.1583214259
-0.7500000000   -0.5000000000   +0.4953780742   -1.1583214259
-0.7500000000   -0.6250000000   +0.4953780742   -0.4922979148
-0.7500000000   -0.6875000000   +0.4953780742   -0.0447964325
-0.7187500000   -0.6875000000   +0.2128474183   -0.0447964325
-0.7031250000   -0.6875000000   +0.0810265592   -0.0447964325
-0.6953125000   -0.6875000000   +0.0173789137   -0.0447964325
-0.6953125000   -0.6914062500   +0.0173789137   -0.0138911236
-0.6933593750   -0.6914062500   +0.0016980959   -0.0138911236
-0.6933593750   -0.6923828125   +0.0016980959   -0.0061079375
[[0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9], [0.30685281944005338, 0.19314718055994662, 0.056852819440053382, 0.068147180559946618, 0.0056471805599466185, 0.025602819440053382, 0.0099778194400533815, 0.0021653194400533815, 0.0017409305599466185, 0.00021219444005338151, 0.00076436805994661849]]

In [7]:
x1 = -1.0
f1 = func(x1)
list_newton = [[0],[abs(x1-x0)]]
print('%-15.10f %+24.25f' % (x1,f1))
for i in range(0, 10):
x1 = x1 - f1 / df(x1)
f1 =func(x1)
print('%-15.10f %+24.25f' % (x1,f1))
list_newton[0].append(i)
list_newton[1].append(abs(x1-x0))

print(list_newton)

-1.0000000000   +3.9049848840251204507012517
-0.7909883534   +0.9068233052059522236731937
-0.7057281263   +0.1025656393923117803979039
-0.6933803632   +0.0018661139743176846650385
-0.6931472621   +0.0000006522702786782019757
-0.6931471806   +0.0000000000000799360577730
-0.6931471806   +0.0000000000000000000000000
-0.6931471806   +0.0000000000000000000000000
-0.6931471806   +0.0000000000000000000000000
-0.6931471806   +0.0000000000000000000000000
-0.6931471806   +0.0000000000000000000000000
[[0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9], [0.30685281944005338, 0.097841172874716609, 0.012580945692322598, 0.00023318267075744803, 8.1533773510500396e-08, 8.659739592076221e-15, 1.3322676295501878e-15, 1.3322676295501878e-15, 1.3322676295501878e-15, 1.3322676295501878e-15, 1.3322676295501878e-15]]

In [8]:
import matplotlib.pyplot as plt

X = list_bisec[0]
Y = list_bisec[1]
plt.plot(X, Y)

X = list_newton[0]
Y = list_newton[1]
plt.plot(X, Y)

plt.yscale("log") # y軸を対数目盛に
plt.show()


## 精度，誤差¶

In [9]:
from decimal import *

def pretty_p(result,a,b,operator):
print('context.prec:{}'.format(getcontext().prec))
print(' %20.14f' % (a))
print( '%1s%20.14f' % (operator, b))
print('-----------')
print( ' %20.14f' % (result))

In [10]:
getcontext().prec = 5

aa = '0.80000'
bb = '3.1415'
cc = '3.1234'
a=Decimal(aa)
b=Decimal(bb)
c=Decimal(cc)
pretty_p(b-c,b,c,'-')

print(b-c)
print(a/(b-c))

context.prec:5
3.14150000000000
-    3.12340000000000
-----------
0.01810000000000
0.0181
44.199

In [11]:
TWOPLACES = Decimal(10) ** -3
getcontext().prec = 4
a=Decimal(aa).quantize(Decimal(10) ** -4)
b=Decimal(bb).quantize(Decimal('0.001'))
c=Decimal(cc).quantize(Decimal('0.001'))

pretty_p(b-c,b,c,'-')

print(b-c)
print(a/(b-c))

context.prec:4
3.14200000000000
-    3.12300000000000
-----------
0.01900000000000
0.019
42.11

In [12]:
TWOPLACES = Decimal(10) ** -2
getcontext().prec = 3
a=Decimal(aa).quantize(Decimal(10) ** -3)
b=Decimal(bb).quantize(Decimal('0.01'))
c=Decimal(cc).quantize(Decimal('0.01'))

pretty_p(b-c,b,c,'-')

print(b-c)
print(a/(b-c))

context.prec:3
3.14000000000000
-    3.12000000000000
-----------
0.02000000000000
0.02
40.0

In [13]:
TWOPLACES = Decimal(10) ** -1
getcontext().prec = 2
a=Decimal(aa).quantize(Decimal(10) ** -2)
b=Decimal(bb).quantize(Decimal('0.1'))
c=Decimal(cc).quantize(Decimal('0.1'))

pretty_p(b-c,b,c,'-')

print(b-c)
print(a/(b-c))

context.prec:2
3.10000000000000
-    3.10000000000000
-----------
0.00000000000000
0.0

---------------------------------------------------------------------------
DivisionByZero                            Traceback (most recent call last)
<ipython-input-13-c01b8733d81c> in <module>
8
9 print(b-c)
---> 10 print(a/(b-c))

DivisionByZero: [<class 'decimal.DivisionByZero'>]

## 最小2乗法¶

In [14]:
import numpy as np
from pprint import pprint

aa = np.array([[-2.0, -561.78952],
[-1.0, -564.14261],
[0.0, -565.47273],
[1.0, -565.8513],
[2.0, -565.36457]])

at = np.transpose(aa)
print(at)

[[  -2.        -1.         0.         1.         2.     ]
[-561.78952 -564.14261 -565.47273 -565.8513  -565.36457]]

In [15]:
%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def f(x, a0, a1, a2):
return a0 + a1*x + a2*x**2

xdata = np.array(at[0])
ydata = np.array(at[1])
plt.plot(xdata,ydata, 'o', color='r')

params, cov = curve_fit(f, xdata, ydata)
print(params)

x =np.linspace(-2,2,20)
y = f(x,params[0],params[1],params[2])
plt.plot(x,y, color='b')

plt.grid()
plt.show()

[ -5.65471459e+02  -8.85879000e-01   4.73656429e-01]


## 常微分方程式¶

In [16]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

def my_plot(xx, vv, tt):
plt.plot(tt, xx, color = 'b', linestyle='--',label="height")
plt.plot(tt, vv, color = 'r', label="velocity")
plt.legend()
plt.xlabel('time')
plt.ylabel('height and velocity')
plt.grid()
plt.show()

def euler2(x0, v0):
v1 = v0 + (-cc * v0- g) * dt
x1 = x0 + v0 * dt
return [x1, v1]

In [17]:
g, dt, cc=9.8, 0.1, 1.0
# tt,xx,vv=[0.0],[0.0],[-10]
tt,xx,vv=[0.0],[1000.0],[0.0]
t = 0.0
for i in range(0,1200):
t += dt
x, v = euler2(xx[-1],vv[-1])
tt.append(t)
xx.append(x)
vv.append(v)

my_plot(xx, vv, tt)

In [37]:
xx[1030],xx[1031]

Out[37]:
(0.3999999999914615, -0.5800000000085379)

In [39]:
vv[1030]

Out[39]:
-9.799999999999994

In [38]:
vv[1030]*3600/1000

Out[38]:
-35.27999999999998
In [ ]: