#!/usr/bin/env python
# coding: utf-8
# # Table of Contents
#
#
#
#
# 2019ペア試験解答例
#
#
#
# cc by Shigeto R. Nishitani 2019
#
# # 微積分
#
# ## 偏微分,saddle point
# 関数$f(x,y)=x^2 - y^2$の2次偏導関数$f_{xx}, f_{xy}, f_{yx}, f_{yy}$を求めよ.
# また,$(x,y)=(0,0)$での判別式
# \begin{equation*}
# D = {f_{xy}(0,0)}^2-f_{xx}(0,0)f_{yy}(0,0) > 0
# \end{equation*}
# を確かめよ.さらに,関数$f(x,y)$をplot3dして鞍点(saddle point)の意味を確認せよ.
# (15点)
#
# In[13]:
from sympy import *
x,y = symbols('x y')
fx = diff(x**2-y**2,x)
fx
# In[18]:
fxx = diff(fx,x)
fxy = diff(fx,y)
print(fxx, ',', fxy)
# In[6]:
fy = diff(x**2-y**2,y)
fy
# In[19]:
fyx = diff(fy,x)
fyy = diff(fy,y)
print(fyx, ',', fyy)
# In[21]:
D = fxy**2 - fxx*fyy
print(D)
# In[11]:
get_ipython().run_line_magic('matplotlib', 'notebook')
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
def f(x,y):
return x**2-y**2 #4*x+2*y-6*x*y
x = np.arange(-3, 3, 0.25)
y = np.arange(-3, 3, 0.25)
X, Y = np.meshgrid(x, y)
Z1 = f(X,Y)
fig = plt.figure()
plot3d = Axes3D(fig)
plot3d.plot_surface(X,Y,Z1)
plt.show()
# Dが正であることが確認できる.plot3dによって3次元表示すると,鞍のような表面となり,極値でないことが確認できる.すなわち,fxxの曲率が正,fyyの曲率は負であり,2方向によって曲率が違い鞍点(saddle point)となる.
# ## (フーリエ積分)
# 関数$f(x)= \sin(x) \sin(2x)$の不定積分を求めよ.$f(x)$を$x=-\pi..\pi$でプロットし,この区間での積分値を求めよ.結果についてコメントせよ. (15点)
#
#
# In[33]:
# kernel->Restart
from sympy import *
init_session()
# In[38]:
integrate(sin(x)*sin(2*x), x)
# In[39]:
get_ipython().run_line_magic('matplotlib', 'inline')
x = Symbol('x')
plot(sin(x)*sin(2*x), (x, -pi, pi))
# In[40]:
integrate(sin(x)*sin(2*x), (x, -pi, pi))
# 不定積分の結果は,ややこしいが,プロットからわかるとおり,関数は正の領域と負の領域が対称に出て来ている.したがって,積分すれば正負で相殺して0となる.
# # 線形代数
# ## ヌルスペース
#
# 行列$A = \left(\begin{array}{cccc}
# 4 & -1 & -1 & 1\\
# 1 & 2 & -1 & -2\\
# \end{array}
# \right)
# $を表現行列とする$R^4\rightarrow R^2$の線形写像$f$のKer($f$)の1組の基底を求めよ.(15点)
#
#
# In[42]:
from sympy import *
init_session()
# In[45]:
A=Matrix([[4,-1,1,1],[1,2,-1,-2]])
A.nullspace()
# ## 対角化
# 行列$A = \left(\begin{array}{cccc}
# 4 & -1 & -1 \\
# 1 & 2 & -1 \\
# 3 & -1 & 0
# \end{array}
# \right)
# $の固有値と固有ベクトルを求めよ.
# また,対角化行列$P$を求めて,
# $P^{-1}AP$と$P^{t}AP$を求め,違いを確かめよ.(15点)
# In[46]:
A = Matrix([[4,-1,-1],[1,2,-1],[3,-1,0]])
A.eigenvects()
# In[47]:
P,D = A.diagonalize()
# In[48]:
P.inv()*A*P
# In[50]:
P.transpose()*A*P
# # 2015年度大学入試センター試験 本試験 数学II・B第2問
#
# 以下のセンター試験問題をpythonでcode化せよ.ただし,関数$f(x)$は
# ``` python
# f =Rational(1,2)*x**2
# f.subs({x:a})
# ```
# などとするべし.
#
# ## 微分係数
# 関数$f(x)=\frac{1}{2}x^2$の$x=a$における微分係数$f'(a)$を
# 求めよう.$h$が0でないとき,
# $x$が$a$から$a+h$まで変化するときの$f(x)$の
# 平均変化率は$\fbox{ ア }+\frac{h}{\fbox{ イ }}$である.
# したがって,求める微分係数は
# \begin{equation*}
# f'(a) = \lim_{h \rightarrow \fbox{ ウ }}
# \left(\fbox{ ア }+\frac{h}{\fbox{ イ }}\right) = \fbox{ エ }
# \end{equation*}
# である.
#
#
#
#
# In[1]:
from sympy import *
init_session()
# In[2]:
a,b,h,x = symbols('a b h x')
f =Rational(1,2)*x**2
f
# In[3]:
f.subs({x:a+h})
# In[4]:
f.subs({x:a+h}) - f.subs({x:a})
# In[5]:
(f.subs({x:a+h}) - f.subs({x:a}))/h
# In[6]:
expand( (f.subs({x:a+h}) - f.subs({x:a}))/h )
# In[7]:
limit(expand( (f.subs({x:a+h}) - f.subs({x:a}))/h ), h, 0)
# 少し冗長ですが,ゆっくりと数式処理のステップを示しています.
#
# * 1/2だと0.5になるので明示的に有理数(Rational)と宣言
# * 関数の代入はsubs(tituion) で実行
# * a+hとaとの差をとる
# * さらに展開(expand)を実行
# * limitを取る.これはテキストでは微分のsectionで紹介している
# ## 放物線,接線
#
# 放物線$y=\frac{1}{2}x^2$を$C$とし,$C$上に点P$\left(a, \frac{1}{2}a^2\right)$をとる.
# ただし,$a > 0$とする.点Pにおける$C$の接線$l$の方程式は
# \begin{equation*}
# y = \fbox{ オ }x - \frac{1}{\fbox{ カ }}a^2
# \end{equation*}
# である.
# 直線$l$と$x$軸との交点$Q$の座標は$\left(\frac{\fbox{ キ }}{\fbox{ ク }}, 0\right)$である.
# 点Qを通り$l$に垂直な直線を$m$とすると,$m$の方程式は
# \begin{equation*}
# y = \frac{\fbox{ ケコ }}{\fbox{ サ }}x + \frac{\fbox{ シ }}{\fbox{ ス }}
# \end{equation*}
# である.
#
#
# ### プロット
# まずは関数をplotしてさらに,aを適当に仮定して図を書いています.
# こいつが面倒なんですが,pythonを許してあげてください.
# コードを覚え必要はなくて,ここにサンプルがあることを覚えておいてください.
# In[8]:
get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
import numpy as np
from sympy import *
def func(x):
return 1/2*x**2
xx = np.linspace(-1, 2.2, 100) #0から2πまでの範囲を100分割した
yy = func(xx)
plt.plot(xx, yy, color = 'b')
plt.plot(2, func(2), "o", color = 'k')
plt.hlines(0, -1, 3, color='k', linestyle='-', linewidth=2)
plt.vlines(0, -1, 3, color='k', linestyle='-', linewidth=2)
plt.grid()
plt.show()
# In[9]:
ll = aa*(x-x0)+y0
# In[10]:
diff(f,x)
# In[11]:
diff(f,x).subs({x:a})
# In[12]:
aa = diff(f,x).subs({x:a})
# In[13]:
x0 = a
y0 = Rational(1,2)*a**2
ll = aa*(x-x0)+y0
expand(ll)
# In[14]:
solve(ll,x)
# In[15]:
mm = bb*(x-x0)+y0
# In[16]:
bb = solve(a*b+1,b)[0]
bb
# In[17]:
mm = bb*(x-a/2)
expand(mm)
# In[18]:
print(mm.subs({a:2}))
ll.subs({a:2})
# In[36]:
get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
import numpy as np
from sympy import *
def func(x):
return 1/2*x**2
def l_func(x):
return 2*x-2
def m_func(x):
return -x/2+1/2
g = plt.subplot()
g.set_ylim([-3,3])
g.set_xlim([-3,3])
g.set_aspect('equal')
xx = np.linspace(-1, 2.2, 100) #0から2πまでの範囲を100分割した
yy = func(xx)
g.plot(xx, yy, color = 'b')
yy = l_func(xx)
g.plot(xx, yy, color = 'r')
yy = m_func(xx)
g.plot(xx, yy, color = 'r')
plt.plot(2, func(2), "o", color = 'k')
plt.text(2+0.5, func(2), "P(a,1/2a^2)", color = 'k')
plt.plot(1,0, "o", color = 'k')
plt.text(1.2, 0.2, "Q(a/2, 0)", color = 'k')
plt.plot(0,1/2, "o", color = 'k')
plt.text(0.2, 1/2+0.2, "A(0, 1/2)", color = 'k')
plt.plot(0,0, "o", color = 'k')
plt.text(0.2, -0.3, "O", color = 'k')
plt.plot(2,0, "o", color = 'k')
plt.text(2.2,-0.3, "R", color = 'k')
plt.hlines(0, -1, 3, color='k', linestyle='-', linewidth=2)
plt.vlines(0, -1, 3, color='k', linestyle='-', linewidth=2)
plt.grid()
plt.show()
# aspect比が1のplotは少しコツがいりそう.上の例を参照してください.
# subplotを使わないとエラーが出たんですが,よくわかりません.
#
# 実はここで間違いをしていました.曲線$m$の傾きを求める式が間違ってたんですが,
# aspect比が1のplotをして,間違いに気がつきました.
# 視覚化はいつも大事ですね.
#
# 直線の方程式ですが,エラーを出しています.
# これはわざとです.直線の方程式の公式
# $$
# y-f(a) = f'(a)(x-a)
# $$
# をpythonでアレンジして
# ``` python
# mm = bb*(x-x0)+y0
# ```
# と書いておくと,どの値が足りないかがわかります.
# 傾き(bb),x0, y0と値(あるいはパラメータ)を
# 指定していって直線が出来上がります.
#
# グラフには以下の問題を解くために,点O,Rを足しています.
# また,後の計算で間違いが起こりにくいように,座標を書き加えています.
#
# ## 面積
# 直線$m$と$y$軸との交点をAとする.三角形APQの面積を$S$とおくと
# \begin{equation*}
# S = \frac{a\left(a^2+\fbox{ セ }\right)}{\fbox{ ソ }}
# \end{equation*}
# となる.また,$y$軸と線分APおよび曲線$C$によって囲まれた図形の面積を$T$とおくと
# \begin{equation*}
# T = \frac{a\left(a^2+\fbox{ タ }\right)}{\fbox{ チツ }}
# \end{equation*}
# となる.
#
# $a>0$の範囲における$S-T$の値について調べよう.
# \begin{equation*}
# S -T = \frac{a\left(a^2-\fbox{ テ }\right)}{\fbox{ トナ }}
# \end{equation*}
# である.
# $a>0$であるから,$S-T>0$となるような$a$ のとり得る値の範囲は
# $ a > \sqrt{\fbox{ ニ }}$である.
# また,$a>0$のときの$S-T$の増減を調べると,
# $S-T$は$a=\fbox{ ヌ }$で最小値$\frac{\fbox{ ネノ }}{\fbox{ ハヒ }}$をとることがわかる.
#
# In[20]:
mm.subs({x:0})
# In[55]:
OAPR = (Rational(1,2)+f.subs({x:a}))*a/2
OAQ = Rational(1,2)*Rational(1,2)*a/2
QPR = (a-a/2)*1/2*a**2/2
OAPR
# In[57]:
SS = factor(OAPR - OAQ - QPR)
SS
# In[61]:
TT = factor(OAPR-integrate(f,(x,0,a)))
TT
# In[64]:
ST = factor(SS-TT)
ST
# In[65]:
plot(ST, (a,0,4))
# In[66]:
solve(ST,a)
# In[67]:
diff(ST,a)
# In[68]:
solve(diff(ST,a),a)
# In[69]:
ST.subs({a:1})
# これでセンター試験の問題は解けました.
# pythonで解いたら綺麗にセンター試験の解答ボックスを
# 埋めていくことができるはずです.
# 出題者は検算をpythonにやらしているのではと訝るほどです.
# でも,次の問題をこのままでは解けません.
# # センター試験数値計算化
#
# 前問の関数を$f(x)=0.49x^2$および放物線$C$の方程式を$y=0.49x^2$として問題を解け.
# 数値解となるので,答えはかっこによらず小数点となる.最後の最小値は-0.08676940ぐらい.(30点)
## 解法の指針
代数計算でできた問題を数値計算にするという作業はよくあります.
研究や実装においては日常的な作業です.
コツは,できるだけ自動化しておくことです.
上の解法では,出て来た結果をそのままベタ打ちしてるところがあります.
そうすると変数が変わった時に「もれ」が出て来ます.
まずは答えを取り出す時に,手打ちではなく変数として取り出す工夫をします.
例えば,P点の座標を$(a,1/2*a^2)$ではなく$(a,f(a))$とするだけで大分自動化が進みます.
他のもやっていきます.
# In[70]:
from sympy import *
init_session()
# In[99]:
a,b,h,x = symbols('a b h x')
#f =Rational(1,2)*x**2
f =0.49*x**2
f
# In[100]:
df = expand( (f.subs({x:a+h}) - f.subs({x:a}))/h )
df
# In[101]:
limit( df, h, 0)
# In[103]:
aa = diff(f,x).subs({x:a})
aa
# In[118]:
x0 = a
y0 = f.subs({x:a})
ll = aa*(x-x0)+y0
expand(ll)
# In[119]:
Q_x = solve(ll,x)[0]
Q_x
# In[120]:
bb = solve(aa*b+1,b)[0]
bb
# In[121]:
mm = bb*(x-Q_x)
expand(mm)
# In[122]:
A_y = mm.subs({x:0})
A_y
# In[123]:
OAPR = (A_y+f.subs({x:a}))*a/2
OAQ = A_y*Q_x/2
QPR = (a-Q_x)*f.subs({x:a})/2
OAPR
# In[124]:
SS = factor(OAPR - OAQ - QPR)
SS
# In[125]:
TT = factor(OAPR-integrate(f,(x,0,a)))
TT
# In[126]:
ST = factor(SS-TT)
ST
# In[127]:
plot(ST, (a,0,4))
# In[128]:
solve(ST,a)
# In[132]:
ST_min = solve(diff(ST,a),a)
ST_min[1]
# In[133]:
ST.subs({a:ST_min[1]})
# どうです.最初は$f=1/2x^2$で結果を確認しながらコードを整理していきます.
# 最後に係数を0.49に変えて実行します.
# 最後の結果が違う人はどこかで変数への変換が不十分なのだと思います.チェックしてみてください.
#
# 「重複」あるいは冗長な部分を省くとみやすくなっていると思います.
# 普通のテキストやネットのサンプルは綺麗なコードを書きたがります.
# でも,ここで示したように,
# 初めから綺麗な答えがあると思わないでください.
# プロもグタグタのcodeを書いて,その後でこそっと直しています.
# 初心者は,特に試行錯誤が必要です.
# この辺りが,「programmingは数学の問題を解くのに似ている」
# と言われる所以(ゆえん)です.
#
# 数式処理のコツがわかったでしょうか?
# プログラミングの極意に似ています.
# 「綺麗な動くコード」です.
# まずは動くコードを書いて,その後重複を省いて綺麗にしてください.
# この過程がREPLと表現されています.
# In[ ]: