2018年の京大数学をプログラミング（Python、Sympy）で解いてみた

こんにちは、現役京大生ブロガーのゲーテ(@goethe_kyodai)です。

普通に解いても面白くないのでプログラミング（Python）で2018年の理系の京大数学を解いてみました。

プログラミングで解くことに重きを置き、厳密な議論は飛ばしていて模範解答とは程遠い解答となっているので高校生のみんなは参考にしないように。

使用言語:Python2系
使用ライブラリ:Sympy,Numpy,Math,matplotlib

Sympyの使い方はこちらを参考にしました。

Sympyの使い方

Pythonの数式処理ライブラリSymPyをWolfram Alpha(Mathematica, Maxima)の代わりに使う方法 - MyEnigma

はてなブログにコードを埋め込むのはこちらを参考にしました。

はてなブログにコードを埋め込む

http://www.tsubasa-note.blog/entry/source-code-highlighting/

問題１
- （１）
- （２）
問題２
問題３
問題４
- 解法１
- 解法２
問題５
- (1)
- (2)
問題６
- (1)
- (2)
感想
参考書

問題１

f:id:a86223990:20180310142205p:plain

（１）

$C_1$ と $C_2$ からyを消去して整理して、

$(b - a)x^2 - 2bx + b + c$ ・・・（１）

を得る。

上式の判別式をdiscriminant()で求めると、

$D = b^2 - (b-a)(b+c)$

となる。

これをsolve()で $b$ について解くと、

$b = \frac{ac}{c-a}$

となる。

これをsubs()で（１）に代入してsolve()でxについて解くと、

$x = \frac{c}{a}$

を得る。

これを $y = ax^2$ にsubs()で代入して、

$t = \frac{c^2}{a}$

を得る。

これらが答えの接点の座標である。

まとめると答えは

$(\frac{c}{a}, \frac{c^2}{a})$

（２）

まず先ほど求めた接点を $(x,y) =(\frac{c}{a}, \frac{c^2}{a})$ と置きます。

そしてsolve()でこれらの連立方程式をx、yについて解くと、

$(a,c) = (\frac{y}{x^2}, \frac{y}{x})$

となります。

これを条件（１）の不等式を等式に置き換えた、 $1 + c^2 -2a$ にsubs()で代入してsolve()でyについて解くと、

$y = -\sqrt(-x^2 + 1) + 1 , \sqrt(-x^2 + 1) + 1$

を得られます。

このyと閉区間 [-1,1] の範囲のxを使って、matplotlibで図示したのが下です。

f:id:a86223990:20180310175749p:plain

本来の条件（１）は不等式なので上図の円に囲まれた部分が答えです。（本来はa,b,cの条件から除外する点があるがめんどくさいので省略）

▶問題１のソースコードを見る

# -*- coding:utf-8 -*-
from sympy import *
from matplotlib import pyplot as plt
import numpy as np
x = Symbol('x')
y = Symbol('y')
a = Symbol('a')
b = Symbol('b')
c = Symbol('c')

#(1)
expr1 = a*x**2 #C1
expr2 = b*(x-1)**2 + c #C2
expr3 = expand(expr2 - expr1) #yを消去

D = Poly(expr3).discriminant() #判別式
b_1 = solve(D,b)[0] #接する条件D=0からbをa,cで表す

expr4 = expr3.subs([(b,b_1)]) #bを消去
contact_x = solve(expr4,x)[0] #xについて解いて
contact_y = a*contact_x**2 #C1にxを代入
print("接点の座標は ({},{})".format(contact_x,contact_y))


#(2)
eq1 = x - contact_x
eq2 = y - contact_y
[a_1,c_1] = solve([eq1,eq2],[a,c])[0] #a,cをy,xのみで表す
expr5 = 1 + c**2 - 2*a_1
expr6 = sympify(expr5.subs([(a,a_1),(c,c_1)]))#a,cを条件(1)に代入
point_y = solve(expr6,y) #yについてとく

#図示
X = np.linspace(-1,1,100)
Y = []

for point_x in X:
    Y.append(point_y[0].subs(x,point_x))
plt.plot(X,Y,'r')
Y = []
for point_x in X:
    Y.append(point_y[1].subs(x,point_x))
plt.plot(X,Y,'r')
plt.xlabel("x-label")
plt.ylabel("y-label")
plt.show()

問題２

f:id:a86223990:20180310142250p:plain

$f(n) = n^3 - 7*n + 9$ と置きます。

因数分解できなさそうなので代入して構造を調べました。

その前に整数nじゃ範囲が広すぎるので、 $f(n)$ が素数になり得るnの範囲を調べておきましょう。

$g(x) = x^3 - 7*x + 9$ と置きます。

これをmatplotlibのpyplotで図示すると下のようになります。

f:id:a86223990:20180310153334p:plain

見た感じ $g(x) > 0$ となるのは $x \geqq -3$ ですね。これで調べる範囲を $n \geqq -3$ に限定できました。

では実際に代入して素数かどうかの判定をとりあえずn = -3 ~ 20 までしてきます。
プログラムで実際に代入した結果がこちら↓

n = -3,f(n) = 3で、「f(n)が素数である」はTrue
n = -2,f(n) = 15で、「f(n)が素数である」はFalse
n = -1,f(n) = 15で、「f(n)が素数である」はFalse
n = 0,f(n) = 9で、「f(n)が素数である」はFalse
n = 1,f(n) = 3で、「f(n)が素数である」はTrue
n = 2,f(n) = 3で、「f(n)が素数である」はTrue
n = 3,f(n) = 15で、「f(n)が素数である」はFalse
n = 4,f(n) = 45で、「f(n)が素数である」はFalse
n = 5,f(n) = 99で、「f(n)が素数である」はFalse
n = 6,f(n) = 183で、「f(n)が素数である」はFalse
n = 7,f(n) = 303で、「f(n)が素数である」はFalse
n = 8,f(n) = 465で、「f(n)が素数である」はFalse
n = 9,f(n) = 675で、「f(n)が素数である」はFalse
n = 10,f(n) = 939で、「f(n)が素数である」はFalse
n = 11,f(n) = 1263で、「f(n)が素数である」はFalse
n = 12,f(n) = 1653で、「f(n)が素数である」はFalse
n = 13,f(n) = 2115で、「f(n)が素数である」はFalse
n = 14,f(n) = 2655で、「f(n)が素数である」はFalse
n = 15,f(n) = 3279で、「f(n)が素数である」はFalse
n = 16,f(n) = 3993で、「f(n)が素数である」はFalse
n = 17,f(n) = 4803で、「f(n)が素数である」はFalse
n = 18,f(n) = 5715で、「f(n)が素数である」はFalse
n = 19,f(n) = 6735で、「f(n)が素数である」はFalse
n = 20,f(n) = 7869で、「f(n)が素数である」はFalse

どうやら $n \geqq 3$ ではf(n)は素数にならないっぽいので、n=3m,3m+1,3m+2(m は整数で１以上)と置いてそれを証明してみます。

とりあえず代入してみます。

f(3m) = $3(9m^3 - 7m + 3)$
f(3m+1) = $3(9m^3 + 9m + 1)$
f(3m+2) = $3(18m^3 + 5m + 1)$

いい感じに3 $\times$ (多項式)になってますね。

あとはf(3m)、f(3m+1)、f(3m+1)を3で割った $9m^3 - 7m + 3$ 、 $9m^3 + 9m + 1$ 、 $18m^3 + 5m + 1$ が２以上の整数になっていれば、f(3m)、f(3m+1)、f(3m+1)が合成数になるので素数にならないことを証明できます。

それをmatplotlibのpyplotで図示して目視で確認しましょう。

f:id:a86223990:20180310161320p:plain

それぞれ２以上になってますね。

以上より答えはn = -3、１、２で解けました。

▶問題２のソースコードを見る

# -*- coding:utf-8 -*-
from matplotlib import pyplot as plt
from sympy import *
import numpy as np
import math

#Prime numberかどうか判定
def is_prime(x):
    if x < 2: return False # 2未満に素数はない
    if x == 2 or x == 3 or x == 5: return True # 2,3,5は素数
    if x % 2 == 0 or x % 3 == 0 or x % 5 == 0: return False # 2,3,5の倍数は合成数
    
    # ためし割り: 疑似素数(2でも3でも5でも割り切れない数字)で次々に割っていく
    prime = 7
    step = 4
    while prime <= math.sqrt(x):
        if x % prime == 0: return False
        
        prime += step
        step = 6 - step
    return True

x = Symbol('x')
n = Symbol('n')
m = Symbol('m')
expr1 = x**3 - 7 * x +9 #与式 f(x)とする

#f(n) > 0 となるnの範囲を図示と目視で調べる
X = np.linspace(-4,4,100)
Y = X**3 - 7 * X + 9
plt.plot(X,Y,label="g(x)")
plt.plot([-4,4],[0,0],label="g(x) = 0")
plt.plot([-3,-3],[-30,40],label="x = -3")
plt.title("g(x) = $x^3 - 7x + 9$")
plt.xlabel("x-LABEL")
plt.ylabel("g(x)-LABEL")
plt.legend()
plt.show()

#n>=-3でf(n)が素数になるかどうか調べる
a = -3
for i in range(24):
    solution = expr1.subs(x,a)
    isPrime = is_prime(solution)
    print("n = {},f(n) = {}で、「f(n)が素数である」は{}".format(a,solution,isPrime))
    a += 1
print("\n")
#f(3m)、f(3m+1)、f(3m+1)の形を調べる
b = [3*m,3*m+1,3*m+2] #m >= 1
for bb in b:
    expr2 = expand(expr1.subs(x,bb))
    print("n = {},f(n) = {}".format(bb,expr2))
print("\n")

#f(3m)/3、f(3m+1)/3、f(3m+1)/3がm>=1で２以上になることを図示して確認
point_x = np.linspace(-3,3,100)
for bb in b:
    expr3 = expand(expr1.subs(x,bb))/3
    print("n = {},f(n)/3 = {}".format(bb,expr3))
    point_y = []
    for point in point_x:
        point_y.append(expr3.subs(m,point))
    plt.plot(point_x,point_y,label = "{}".format(expr3))

plt.plot([-3,3],[2,2],label = "f/3 = 2")
plt.plot([1,1],[-400,400],label = "m = 1")
plt.xlabel("m-label")
plt.title("$f(3m)/3,f(3m+1)/3,f(3m+2)/3$")
plt.legend()
plt.show()

問題３

f:id:a86223990:20180310142256p:plain

f:id:a86223990:20180310180621p:plain

変数 $\theta (0 < \theta < \alpha)$ を導入して、条件から上図のように角を表します。

外接円の半径が１であることと正弦定理より
$\frac{AB}{sin(\pi - ( \alpha + \theta)) } = \frac{BC}{sin(\theta)}$
$= \frac{CD}{sin(\alpha - \theta)} = \frac{DA}{sin(\theta)} =2$

よって
AB = $2sin(\pi - ( \alpha + \theta))$
BC = $2sin(\theta)$
CD = $2sin(\alpha - \theta)$
DA = $2sin(\theta)$
と表せます。

これらを使ってkは

k = $16sin(\theta)^2sin(\pi - ( \alpha + \theta)sin(\alpha - \theta)$

と表せます。（ここまでプログラミング登場なし）

入試問題的には最大値を $\alpha$ で表せばいいのですが、プログラミングを駆使して解きたいので0 < $\alpha \leqq \frac{\pi}{2}$ を満たす任意の実数の入力に対して、線形探索で最大値を出すプログラムを書きました。

そのプログラムのアルゴリズムは
１： $\alpha$ を受け取る
２：0 < $\theta$ < $\alpha$ の範囲で任意個に均等に分けた $\theta$ を用意する。
３：さっき分けた $\theta$ を小さい方から順に $\alpha$ と共にkに代入して計算し、kの最大値を線形探索で求める

線形探索なので計算量はO(n)(n は $\theta$ を分割した数)です。

▶問題３のソースコードを見る

# -*- coding:utf-8 -*-
from sympy import *
import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

a = Symbol('a') #めんどいのでαをaで表現
b = Symbol('b') #めんどいのでthetaをbで表現

#立式
AB = 2 * sin(pi - (a+b))
BC = 2 * sin(b)
CD = 2 * sin(a-b)
DA = 2 * sin(b)
k = expand(AB*BC*CD*DA)

#aの値に応じてkの最大値を求める
constant_a = math.pi/4.0 #ここに0 < constant_a <= pi/2を満たす好きなA(alpha)を入力
roop = 10000 #0 < B < Aの範囲でBをroop個に均等に分ける
            #roopが大きいほど答えの精度がよくなる
B = np.linspace(0,constant_a,roop+2)
B = np.delete(B,0,axis=None)    #0 < theta < alphaより端点(0,alpha)を除去
B = np.delete(B,-1,axis=None)
max_k = -16     #alpha = Aの時のkの最大値
max_b = B[0]    #kがmaxを取る時のb
k = k.subs([(a,constant_a)])
for bb in B:
    tmp = k.subs([(b,bb)])
    if tmp > max_k:
        max_b = bb
        max_k = tmp
print("alpha = {},theta = {}の時、kの最大値は{}".format(constant_a,max_b,max_k))

問題４

f:id:a86223990:20180310142303p:plain

解法１

$z_n = 1$ となる確率をP( $z_n$ = 1)とおく。

P( $z_n$ = 1 ) = $\frac{z_n = 1となるz_nの数}{全ての z_n}$ 　です。

[全ての $z_n$ ]は $2^n$ で、[ $z_n$ = 1となる $z_n$ の数] ]は再帰を使った全探索で求められます。

計算量がO( $2^n$ )になってしまいますが、計算時間の制限がなければ任意の入力nに対して一応この解法で解けます。

▶問題４解法１のソースコードを見る

# -*- coding:utf-8 -*-
from sympy import *
import math

expr1 = -0.5 + sqrt(3)*I/2
n = 10 #ここに好きなnを入れる
count = 0 #z(n) == 1となるz(n)の数

#depth番目の数列zを再帰で求める
def calcz(z,depth):
    global n
    global count
    
    if depth == n:
        if z == 1:
            count +=1
    else:
        #表の場合
        front = expand(expr1*z)
        calcz(front,depth+1)
        
        #裏の場合
        back = expand(conjugate(z))
        calcz(back,depth+1)

#初期値z1 = expr1,1でそれぞれzn求める
calcz(expr1,1)
calcz(1,1)

prob = (count + 0.0)/2**n
print("n = {},probability = {}".format(n,prob))

解法２

解法１は計算量が大きすぎるので速い方法で解いてみます。

a = 1、b = $\frac{-1+\sqrt{3} i}{2}$ 、c = $\frac{-1-\sqrt{3} i}{2}$ と置きます。

$z_n$ がa、b、cとなる確率を $a_n$ 、 $b_n$ 、 $c_n$ と置きます。

それらを用いた以下の漸化式が成り立ちます。

$a_n$ + $c_n$ = 2 $a_{n+1}$
$a_n$ + $c_n$ = 2 $b_{n+1}$
$b_n$ = $c_{n+1}$

A = $\left( \begin{array}{ccc} 1 & 0 & 1 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{array} \right)$ 、B = $\left( \begin{array}{ccc} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 1 \end{array} \right)$ 、 $x_{n} = (a_n,b_n,c_n)$ とすると上の漸化式は

A $x_n$ = B $x_{n+1}$

と表せます。

式変形すると

$x_{n}$ = $B^{-1}Ax_{n-1}$

と表せます。C = $B^{-1}A$ と置くと

$x_{n}$ = C $x_{n-1}$

となるので

$x_{n}$ = $C^{n-1} x_1$

と表せます。

$C^{n-1}$ をプログラミングの計算力を使ってCをn-1回掛けて求めてもいいですが、計算量がO( $(3^2)^n$ )になってしまうので対角化しましょう。

D = $P^{-1} C P$ となるようにCを対角化します。この時、Cの固有値をdiagnolize()で求めて、 $\lambda_1、\lambda_2、\lambda_3$ とすると、

$C^{n-1}$ = $P D^{n-1} P^{-1}$ = $P \left( \begin{array}{ccc} \lambda_1^{n-1} & 0 & 0 \\ 0 & \lambda_2^{n-1} & 0\\ 0 & 0 & \lambda_3^{n-1} \end{array} \right) P^{-1}$

と表せて計算量がO(n)ですみます。

この $C^{n-1}$ と $x_1 = (\frac{1}{2},\frac{1}{2},0)$ を使って $x_n = (a_n,b_n,c_n)$ を求められます。

この $a_n$ が求める確率P( $z_n$ = 1)です。

▶問題４解法２のソースコードを見る

# -*- coding:utf-8 -*-
from sympy import *
import math
n = Symbol('n')

A = Matrix([[1,0,1],[1,0,1],[0,1,0]])
B = Matrix([[2,0,0],[0,2,0],[0,0,1]])
C = B.inv()*A
(P,D) = C.diagonalize() #対角化 D = P^-1CP

#Dのn-1乗のDDを作る
DD = zeros(3) #Dのn-1乗
i = 0
while(i <= 3*3 - 1):
    DD[i] = D[i]**n
    i += 3 + 1

#x_nをx_n = P(Dのn-1乗)P^-1*x_1で計算
x_1 = Matrix([0.5,0.5,0]) #x_1 = (a_1,b_1,c_1)
x_n = P*DD*P.inv()* x_1 #x_n = (a_n,b_n,c_n)
print("a_n= {}".format(x_n[0]))

問題５

f:id:a86223990:20180310142307p:plain

(1)

$u(t)、v(t)$ を普通に求めると、

$(u(t), v(t))$

$= (t + \frac{1}{\sqrt{t^2 + 1}},-\frac{t}{\sqrt{t^2 + 1}} + log(t))$

これらをdiff()を使ってtで微分すると答えを出せます。

答えは

$(\frac{du}{dt}, \frac{dv}{dt})$
$= (1 - \frac{t}{(t^2+1)\sqrt{t^2 + 1}},\frac{1}{t} - \frac{1}{(t^2+1)\sqrt{t^2 + 1}})$

です。

(2)

$L_1(r) = \int_1^r \sqrt{1+\frac{1}{t^2}} dt$ 、 $L_2(r) = \int_1^r \sqrt{ ( \frac{du}{dt})^2+( \frac{dv}{dt} )^2 } dt$ と表せます。

integrate()で積分して $L_1(r)$ と $L_2(r)$ を求めてlimit()で $\lim_{ r \to \infty} L_1(r) - L_2(r)$ を求めようと試みましたが、、、

$L_1(r)$ はすぐ求められました、 $L_2(r)$ の式がゴツすぎるせいか中々結果を吐き出してくれませんでした。

なので式を簡単化します。

(1)で求めた導関数をよく見ると、 $\frac{du}{dt}$ = $t \times \frac{dv}{dt}$ の関係になってます。

これを使うと $L_2(r) = \int_1^r \frac{du}{dt} \sqrt{ 1+( \frac{1}{t^2} ) } dt$ と表せ、さらにこれを使って $L_1(r) - L_2(r)$ を簡単化して

$L_1(r) - L_2(r) = \int_1^r \frac{1}{1+t^2} dt$

と表せます。（Sympyのcancel()やsimplify()を使って $\frac{du}{dt}$ = $t \times \frac{dv}{dt}$ からここまで簡単化できると思ったんですが中々できなくて手で計算しました。）

このように式変形してlim()でやっと解けるようになります。

▶問題5のソースコードを見る

# -*- coding:utf-8 -*-
from sympy import *
import math
t = Symbol('t')
r = Symbol('r')

#u,v立式
u = t + 1/sqrt(t**2 + 1)
v = log(t) - t/sqrt(t**2 + 1)
print("(u(t),v(t)) = ({}, {})\n".format(u,v))

#u,vをtで微分
du = diff(u,t)
dv = diff(v,t)
print("(du/dt,dv/dt) = ({}, {})\n".format(du,dv))

expr1 = 1/(t**2+1)
expr2 = integrate(expr1,t).subs(t,1) - integrate(expr1,t).subs(t,r)
solution = limit(expr2,r,0)
print("lim(L1(r) - L2(r)) = {}".format(solution))

問題６

f:id:a86223990:20180310142313p:plain

(1)

まずD= (0,0,0)、A=( $a_1,a_2,a_3$ )、B=( $b_1,b_2,b_3$ )、C= ( $c_1,c_2,c_3$ )と置きます。

その座標を使って計算したPQ•ABをLと置きます。PQとABが垂直になることを示せばいいのでL=0を示せばいいですね。

条件AC = BD、AD = BCをsolve()で解いて $a_1、a_2をa_3、b_1、b_2、b_3、c_1、c_2、c_3$ で表して、それらをLにsubs()してsymplify()とexpand()してもLが0にならなかったので方針を変更しました。

条件AC = BD、AD = BCを二乗した $AC^2 = BD^2、AD^2 = BC^2$ を辺々足してexpand()で展開すると、

$2a_1 ^2 - 2a_1 b_1 -2a_1 c_1 + 2a_2 ^2 - 2a_2 b_2$

$- 2a_2 c_2 + 2a_3 ^2 -2a_3 b_3 - 2a_3 c_3 = 0$

になります。左辺をMと置きます。

Lをexpand()で展開すると、

$L$

$= -0.5a_1 ^2 + 0.5a_1 b_1 +0.5a_1 c_1 - 0.5a_2 ^2$

$+ 0.5a_2 b_2 + 0.5a_2 c_2 - 0.5a_3 ^2 +0.5a_3 b_3 + 0.5a_3 c_3$

となります。

二つの式を目視で比べると、L = - 0.25 Mとなっていることがわかるので、L = 0となり、証明できました。(結局目視)

▶問題6(1)のソースコードを見る

# -*- coding:utf-8 -*-
from sympy import *
a1 = Symbol('a1')
a2 = Symbol('a2')
a3 = Symbol('a3')
b1 = Symbol('b1')
b2 = Symbol('b2')
b3 = Symbol('b3')
c1 = Symbol('c1')
c2 = Symbol('c2')
c3 = Symbol('c3')

BA = Matrix([a1,a2,a3])
BD = Matrix([b1,b2,b3])
BC = Matrix([c1,c2,c3])
AC = BC - BA
AD = BD - BA
PA = 0.5*BA
AQ = 0.5*AD + 0.5*AC
PQ = PA + AQ

#条件から
expr1 = AC.transpose()*AC - BD.transpose()*BD #a,b,c
expr2 = AD.transpose()*AD - BC.transpose()*BC #a,b,b
M = expr1 + expr2
print("M = {}".format(M[0]))
L = PQ.transpose()*BA
print("L = {}".format(L[0]))