はじめに
物理学科の2年生にPythonを教えています。最初に断っておきますがもし、この記事で扱うプログラミングのレベルがとても低いと感じたらそれは学生が悪いんじゃなく僕の教育スキルが低いんだと思ってください。記事の主旨はあくまで生成AI以前と以後で同じ試験ができなくなってどういう工夫をしたかというノウハウ（？）の共有です。
生成AI以前
生成AIが出てくる前はこんな試験でした。
ヤコビ法で行列の対角化をして、対角化した行列と対角化に必要な直交行列を出力するブログラムをPython で作成せよ。対角化する行列は当日の試験時間内に指示する
これを持ち込み制限なし、インターネットアクセスも許すという環境で課していました（試験時間60分）。勿論、ヤコビ法のライブラリとかの使用は禁止で、自分の能力でスクラッチから書くことが条件です。簡単すぎると思うかもしれませんが、物理学科の学生は主に物理の勉強に時間を割いていることが多く、プログラミングにはそれほど時間を割いて勉強しないのが普通です。なので講義で3×3の行列の例をやったとしてもこれが4×4に成っちゃうと間違っちゃうみたいなレベルの人もたくさんいます。ここに評価がありますが、充実度4/5、楽単度4/5なのでまあまあ学生には評価されてたってことでしょう（知らんけど）。
しかし、この問題をChatGPT(GPT-4)に食わすと
ヤコビ法を用いて行列を対角化するPythonプログラムを作成します。この方法は、与えられた対称行列を直交変換によって対角化し、結果として対角化された行列とその変換に使用される直交行列を出力します。
以下はそのプログラムの例です。まず、適当な対称行列を用意し、その行列をヤコビ法で対角化する関数を定義します。その後、対角化された行列と直交行列を出力します。
このプログラムでは、対称行列が必要です。もし具体的な行列があれば、それを使用することができます。それでは、プログラムを作成しましょう。
ヤコビ法を用いて行列を対角化するプログラムを作成し、実行しました。以下は、例として与えられた対称行列とその対角化の結果です。
対称行列 $A$ は次のようになります：
A = \begin{bmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{bmatrix}
対角化された行列 $D$ と直交行列 $V$ は次のようになります：
D = \begin{bmatrix} 2.00001579 & 0 & 0 \\ 0 & 2.60876925 & 0 \\ 0 & 0 & 1.39121496 \end{bmatrix}
V = \begin{bmatrix} 0.69824165 & -0.69828068 & 0.15767904 \\ 0.00044808 & 0.2206908 & 0.97534372 \\ -0.71586199 & -0.68095496 & 0.1544084 \end{bmatrix}
ここで、$D$ は元の行列 $A$ の対角化された形であり、$V$ はそれに対応する直交行列です。この直交行列を用いて、$A$ を対角化することができます。
import numpy as np
def jacobi_eigensystem(A, tol=1e-10, max_iterations=100):
"""
Perform Jacobi rotation algorithm to diagonalize a symmetric matrix.
Args:
- A: The symmetric matrix to be diagonalized
- tol: Tolerance for the off-diagonal elements
- max_iterations: Maximum number of iterations
Returns:
- D: Diagonalized matrix
- V: Matrix of eigenvectors (orthogonal matrix)
"""
n = A.shape[0]
V = np.eye(n)
D = np.copy(A)
for _ in range(max_iterations):
# Find largest off-diagonal absolute element
i, j = np.unravel_index(np.argmax(np.abs(D - np.diag(np.diagonal(D)))), D.shape)
if np.abs(D[i, j]) < tol:
break
# Calculate the rotation angle
theta = np.arctan2(2 * D[i, j], D[j, j] - D[i, i]) / 2
# Create the rotation matrix
P = np.eye(n)
P[i, i] = np.cos(theta)
P[j, j] = np.cos(theta)
P[i, j] = -np.sin(theta)
P[j, i] = np.sin(theta)
# Update D and V
D = P.T @ D @ P
V = V @ P
return np.diag(np.diagonal(D)), V
# Example symmetric matrix
A = np.array([[2, -1, 0],
[-1, 2, -1],
[0, -1, 2]])
# Diagonalize A using Jacobi method
D, V = jacobi_eigensystem(A)
D, V
みたいな有様でさすがにこれでは試験にならないよなあ、と思ったわけです。
今年の問題
こんな感じにしてみました。
関数$f(x)$を0からyまでニュートン・コーツの公式（近似に用いるn次多項式のnは２以上とする）で積分した関数を$F(y)$とする。
\frac{dz}{dy} = z*F(y)
という微分方程式をルンゲ・クッタの公式で積分するプログラムをPythonで作れ。
但し$z(0)$の値は勝手に決めてよい。
但し、関数$f(x)$には具体的な関数を入れて計算結果が正しいことを必ず確認せよ。
但しニュートン・コーツ法、ルンゲ・クッタ法を実行するライブラリなどは使わず、自分でプログラムを書くこと。
ヒント：
\frac{dz}{dy}=z*F(y)
は以下の様にして簡単に変数分離法で解くことができる
\int \frac{dz}{z} = \log(z)= \int F(y) dy
z = C*\exp \left \lbrace \int F(y) dy \right \rbrace
$C$は$z(0)$によって決まる定数。したがってF(y)が積分できるようにf(x)を決めておけば正解は容易に手計算できるので結果が一致することを確認するのは難しくない。
ニュートンコーツとルンゲクッタは講義で説明してあり、その時に紹介した模範解答も自由にアクセスできる状態で試験を受けます。その条件ならChatGPTを使わなくても簡単だろみたいな話は置いておきます。この記事の目的はあくまでChatGPTでは正解できない問題をあんまり問題そのものを難しくせずに作るにはどうするか、ということだからです。
生成AIの結果
プロンプトを工夫すればこの問題も解けるんでしょうけどあまり工夫せずに素で入力してみます。
関数$f(x)=x**2$を0からyまでニュートン・コーツの公式（近似に用いるn次多項式のnは２とする）で積分区間の分割数は100で積分した関数を$F(y)$とする。
$dz/dy = z*F(y)$
という微分方程式をルンゲ・クッタの公式で積分するプログラムをPythonで作れ。
但しz(0)=1とする。
但しニュートン・コーツ法、ルンゲ・クッタ法を実行するライブラリなどは使わず、自分でプログラムを書くこと。
出力はこんな感じ。
了解しました。この問題を解決するために、以下の手順を実行します。
- ニュートン・コーツの公式を用いて、関数 $ f(x) = x^2 $ を区間 $[0, y]$ で近似的に積分します。この場合、2次の多項式（$n = 2$）を使用し、積分区間は100の分割数で扱います。これによって得られる積分値を $ F(y) $ とします。
- 微分方程式 $ \frac{dz}{dy} = z \cdot F(y) $ をルンゲ・クッタ法（4次）を用いて数値的に解きます。初期条件は $ z(0) = 1 $ です。
まずは、ニュートン・コーツの公式を用いて $ F(y) $ を計算する関数を実装しましょう。次に、この $ F(y) $ を使用してルンゲ・クッタ法で微分方程式を解く関数を実装します。それでは、プログラムを書き始めましょう。
プログラムを実行した結果、関数 $ f(x) = x^2 $ をニュートン・コーツの公式で近似的に積分し、その結果を用いて微分方程式 $ \frac{dz}{dy} = z \cdot F(y) $ をルンゲ・クッタ法（4次）で数値的に解きました。この計算により、様々な $ y $ の値に対する $ z $ の値を得ました。
結果は $ y $ の値とそれに対応する $ z $ の値のペアとして出力されています。この結果を使って、必要に応じて更なる分析やプロットを行うことができます。ご希望があれば、これらの値をグラフにプロットすることも可能です。ご要望をお知らせください。
import numpy as np
def newton_cotes_integration(y, n=100):
"""
Perform Newton-Cotes integration of the function f(x) = x**2 over the interval [0, y].
The interval is divided into 'n' equal parts.
"""
# Define the function to integrate
f = lambda x: x**2
# Create equally spaced points in the interval [0, y]
x = np.linspace(0, y, n + 1)
# Compute the width of each subinterval
h = (y - 0) / n
# Compute the integral using the Newton-Cotes formula (Trapezoidal rule)
integral = h * (0.5 * f(x[0]) + np.sum(f(x[1:-1])) + 0.5 * f(x[-1]))
return integral
def runge_kutta_integration(f, y0, z0, h, steps):
"""
Solve the differential equation dz/dy = f(y, z) using the 4th order Runge-Kutta method.
Initial condition: z(y0) = z0.
"""
# Initialize array for y and z values
y_values = np.linspace(y0, y0 + steps * h, steps + 1)
z_values = np.zeros(steps + 1)
z_values[0] = z0
# Runge-Kutta method
for i in range(steps):
k1 = h * f(y_values[i], z_values[i])
k2 = h * f(y_values[i] + 0.5 * h, z_values[i] + 0.5 * k1)
k3 = h * f(y_values[i] + 0.5 * h, z_values[i] + 0.5 * k2)
k4 = h * f(y_values[i] + h, z_values[i] + k3)
z_values[i + 1] = z_values[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6
return y_values, z_values
# Define the differential equation dz/dy = z * F(y)
def differential_eq(y, z):
F_y = newton_cotes_integration(y)
return z * F_y
# Initial conditions
y0 = 0
z0 = 1
h = 0.1 # Step size
steps = 100 # Number of steps
# Solve the differential equation using Runge-Kutta method
y_values, z_values = runge_kutta_integration(differential_eq, y0, z0, h, steps)
y_values, z_values # Return the computed values for further analysis or plotting
一見よさそうですが、実は積分が2次のニュートンコーツじゃなくただの台形公式になってしまっています。
# Compute the integral using the Newton-Cotes formula (Trapezoidal rule)
integral = h * (0.5 * f(x[0]) + np.sum(f(x[1:-1])) + 0.5 * f(x[-1]))
これに気付かずに提出した学生は見事に０点になるので生成AIでそのまま作ってコピペしても点が取れない問題を作ることができました!
ちなみにChatGPT(GPT-4)が2次のニュートンコーツができないわけじゃなく
関数$f(x)=x**2$を0からyまでニュートン・コーツの公式（近似に用いるn次多項式のnは２とする）で積分区間の分割数は100で積分した関数を$F(y)$とする。
このプログラムをpythonでライブラリを使わずに作れ
と入れると
The function $F(y)$, which represents the integral of $f(x) = x^2$ from 0 to $y$ using the Simpson's rule (a Newton-Cotes formula with $n=2$), with the integration interval divided into 100 subintervals, has been successfully implemented in Python.
For example, if we evaluate $F(y)$ at $y = 5$, the result is approximately 41.67.
import numpy as np
def simpson_integration(func, a, b, n):
"""
Apply Simpson's rule for numerical integration.
:param func: Function to integrate.
:param a: Lower limit of integration.
:param b: Upper limit of integration.
:param n: Number of subintervals (must be even).
:return: Approximated integral of the function.
"""
if n % 2 == 1:
raise ValueError("Number of subintervals n must be even.")
h = (b - a) / n
x = np.linspace(a, b, n + 1)
y = func(x)
integral = h/3 * np.sum(y[0:-1:2] + 4*y[1::2] + y[2::2])
return integral
def f(x):
return x ** 2
def F(y):
return simpson_integration(f, 0, y, 100)
# Test the function for a specific value of y, for example y=5
y_test = 5
result = F(y_test)
result
みたいな感じでちゃんと2次のニュートンコーツを作ることができます（答えが英語になってしまいましたが）。問題文をちょっと複雑にしたことでChatGPTではうまく対応できなくさせることができ、コピペで回答したのに満点になることを避けることができるようになりました。
おわりに
これがいい方法かどうかは解りませんが、問題をちょっとひねってやるとコピペでは解けない問題を作れるみたいです。皆さんの工夫も教えてもらえると嬉しいです。
ただ、例年に比べて明らかに問題が難しくなっているので学生が解けるかどうかですね。ただ、まあ、ChatGPTと同じことしかできない学生にはプログラミングの単位をあげられなくなったと考えれば合理的なのかもしれません。ただ、ChatGPTが出る前なら単位があったであろう学生が取れなくなったと思うと気の毒な気もするのですが...。
最後にしつこい様ですが、物理学科にも普通にプログラミングに卓越した学生はたくさんいますのでこれをみて弊学科の学生のプログラミング能力の標準であると思われてしまうと彼らに失礼なのでよろしくお願いします。