#!/usr/bin/env python
# coding: utf-8

# <h1>Table of Contents<span class="tocSkip"></span></h1>
# <div class="toc"><ul class="toc-item"><li><span><a href="#試験解答例18" data-toc-modified-id="試験解答例18-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>試験解答例18</a></span><ul class="toc-item"><li><span><a href="#簡単な行列計算" data-toc-modified-id="簡単な行列計算-1.1"><span class="toc-item-num">1.1&nbsp;&nbsp;</span>簡単な行列計算</a></span></li><li><span><a href="#数値解の収束性" data-toc-modified-id="数値解の収束性-1.2"><span class="toc-item-num">1.2&nbsp;&nbsp;</span>数値解の収束性</a></span></li><li><span><a href="#精度，誤差" data-toc-modified-id="精度，誤差-1.3"><span class="toc-item-num">1.3&nbsp;&nbsp;</span>精度，誤差</a></span></li><li><span><a href="#最小2乗法" data-toc-modified-id="最小2乗法-1.4"><span class="toc-item-num">1.4&nbsp;&nbsp;</span>最小2乗法</a></span></li><li><span><a href="#常微分方程式" data-toc-modified-id="常微分方程式-1.5"><span class="toc-item-num">1.5&nbsp;&nbsp;</span>常微分方程式</a></span></li></ul></li></ul></div>

# # 試験解答例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)


# In[2]:


import scipy.linalg as linalg
inv_a = linalg.inv(aa)
pprint(inv_a)


# In[3]:


pprint(np.dot(inv_a,aa))


# ## 数値解の収束性

# 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)


# 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()


# 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)


# 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)


# 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))


# 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))


# 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))


# 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))


# ## 最小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)


# In[15]:


get_ipython().run_line_magic('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()


# ## 常微分方程式

# In[16]:


get_ipython().run_line_magic('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]


# 地面につくのは，1030 x 0.1 = 103秒後

# In[39]:


vv[1030]


# 速度は，-9.8 m/sec，時速に直すと35.3km/hour．

# In[38]:


vv[1030]*3600/1000


# In[ ]: