生成AI時代のプログラミング科目の試験問題の試み（その２）ではo1でも解けない期末問題をやってみましたという話をしました。
話題のDeepSeek-r1でやってみました。この記事にもあるようにDeepSeek-V3で「DeepThink」を指定するだけでDeepSeek-r1は使えます。結果は以下の通り。

思考過程

はい、この問題を解決するために、まずは全体の流れを整理してみます。問題は、微分方程式 dy/dx = f(x, y) をルンゲクッタ法で解き、得られた解を10×10の行列Aにし、さらにベクトルbとの連立方程式Ax=bをガウスの消去法で解くプログラムを作成することです。ただし、Numpyは使っても良いですが、既存のルンゲクッタやガウス消去法のパッケージは使わないようにします。

まず、ステップごとに考えていきましょう。最初に、微分方程式をルンゲクッタ法で解く部分を作成する必要があります。次に、その解を10×10の行列Aに変換します。そして、ガウスの消去法を実装して、連立方程式を解きます。

まず、ルンゲクッタ法の実装について考えます。ルンゲクッタ法は、4次の公式を使うことが一般的ですね。ステップ幅hを決めて、x=1からx=100まで解を求めます。ここで、yがベクトルである場合を考慮する必要があります。問題では、Aが10x10行列なので、yはおそらく10次元のベクトルでしょう。そのため、微分方程式の右辺f(x, y)は、yを10次元ベクトルとして受け取り、同じ次元のdy/dxを返す関数にする必要があります。ただし、問題文では具体的なf(x, y)や初期値y(0)は任意に決めてよいとありますので、適当な関数を定義します。例えば、f(x, y) = -y のような単純な例を考えるかもしれませんが、その場合の解は指数関数になります。しかし、行列Aを作る部分がどのようになるのか、もう少し考える必要があります。

ここで疑問が湧きます。問題文ではx=1からx=100までの解yを10×10の行列にする、とありますが、これはどういう意味でしょうか？おそらく、xの範囲を100分割して、各ステップのyの値を取るのではなく、x=1からx=100までの各x点（例えば、x=1, 2, ..., 100）に対応するyの値を10×10の行列に配置する、という解釈かもしれません。でも100個のx点があるのに、行列が10×10なので100要素必要です。つまり、xを1から100まで100個の点で計算し、各点でのy（10次元ベクトル）を縦に並べて、行列Aを作るということでしょうか？しかし、行列Aのサイズが10x10なので、各x点でのyの値を何らかの形で10x10にまとめる必要があります。この部分の解釈が少し曖昧です。

例えば、解を求める際に、x=1からx=100までを解く過程で、ステップ数をどう設定するか。例えば、ステップ幅h=1として、x=1,2,...,100まで100ステップ進む場合、各ステップのyの値が10次元ベクトルだとすると、その100個のベクトルをどのように10x10行列に配置するのか。しかし、100個のベクトルがあると、10x10行列に入れるには各ベクトルを1行または1列に配置する必要がありますが、それでは10x10行列に収まらない。ここで矛盾が生じます。

おそらく、問題の設定に誤解があるかもしれません。もしくは、xの範囲が1から100で、その間を離散化して10点だけ取るのでしょうか？例えば、xのステップを等間隔で10分割して、各点でのyの値を10次元ベクトルとして、それらを並べて10x10行列を作る。そうすると、x=1からx=100までを10ステップで進む場合、ステップ幅は(100-1)/9 = 11となりますが、それは現実的ではないかもしれません。あるいは、微分方程式の解を求める際に、各xの値（例えばx=1,2,...,10）に対してyの値を求め、その各yの値が10次元ベクトルなら、それらを10x10の行列として組み立てるという可能性もあります。たとえば、x=1からx=10までの各点で、yが10次元ベクトルであれば、各xに対応するyを横に並べると10x10の行列ができます。この場合、x=1からx=10までの解を求め、各xでのyを列として行列Aを作る。そして、Ax=bを解く。という流れになるかもしれません。しかし、問題文ではx=1からx=100までの解を10x10行列にすると書かれています。100個のx点をどう10x10=100要素にするか、という点では、xを1から100まで100個の点で計算し、各点のy（スカラー？）を10x10行列に入れる、という可能性もありますが、そうするとyは各xに対してスカラー値でなければなりません。しかし、Ax=bは連立一次方程式なので、Aが10x10行列で、xとbが10次元ベクトルです。つまり、問題の設定がやや不明確です。

ここで、問題文を再確認します。「x=1からx=100までの解yを１０×１０の行列にしたものをAとする」とあります。この「解y」が各xに対応する値であるなら、例えばxが1から100までの各点でyがスカラー値の場合、100個のyを10x10の行列に並べる必要があります。しかし、そうすると、行列Aは10x10で、ベクトルbは長さ10なので、Ax=bの方程式を解く際に、Aが10x10、xが10次元、bが10次元となり、これは整合性があります。しかし、もしyが各ステップで10次元ベクトルであれば、行列Aは各ステップのyを並べると、100ステップの場合、100次元のベクトルになりますが、これは10x10の行列に収まらないので矛盾します。したがって、おそらく、yは各x点でスカラー値をとり、x=1から100までの100個のyの値を10x10の行列Aに格納するという解釈が正しいのではないかと考えられます。ただし、この場合、Ax=bはAが10x10行列、xが10次元ベクトル、bが10次元ベクトルなので、数学的に可能です。しかし、元の微分方程式はdy/dx = f(x,y)であり、yがスカラーである場合です。この解釈で進めます。

ただし、別の可能性として、微分方程式の解yが各ステップで10次元ベクトルであり、x=1からx=100までの各x点（100個）でそのベクトルを取得し、それらを10x10行列に配置する方法が考えられます。しかし、100個のベクトルを10x10行列に配置するには、各ベクトルが10要素ある場合、例えば、10個のx点を選んで、それぞれの点での10次元ベクトルを列として並べることで10x10行列を作るという方法があります。つまり、x=1からx=100までを10個の等間隔の点でサンプリングし、各点のyベクトルを列として並べます。例えば、x=1, 11, 21,...,91 の10点を選び、各点でのy（10次元）を列として並べると、10x10の行列Aができます。この解釈であれば、プログラム的には、微分方程式を解く際にステップ数を調整し、10個のサンプル点を選ぶ必要があります。しかし、問題文の記述が曖昧で、どのように解釈すべきか迷います。

ここで、問題文を可能な限り解釈して進める必要があります。まず、微分方程式を解く部分。次に、その解を10x10行列Aに変換。そして、ガウスの消去法でAx=bを解く。

まず、微分方程式の解き方。例として、yの初期値をy(1) = [1, 0, ..., 0]^Tのような10次元ベクトルと仮定し（x=1が初期値？問題文ではx=1からx=100までの解を求めるとあるので、初期条件はx=1のときのyの値が与えられていると考えられます）、微分方程式dy/dx = f(x, y)を解いて、x=100まで進めます。そして、その解を10x10の行列Aに格納する必要があります。例えば、x=1からx=100までの間を等間隔で10分割し、各ステップのyの値を記録して行列を作成します。ステップ数が10ステップなら、x=1, 11, 21,...,100（合計10点）となり、各点でのyの値（10次元ベクトル）を列として並べると10x10行列Aができます。この場合、Aの各列が各x点でのyベクトルとなります。そうすると、Ax=bは、Aが10x10行列、xが10次元ベクトル、bも10次元ベクトルとなるので、問題なく解けます。

では、具体的にプログラムの流れを考えます。

微分方程式dy/dx = f(x, y)をルンゲクッタ法で解く。ここでyは10次元ベクトル。
解を求める範囲はx=1からx=100までで、この間を10ステップ（間隔h= (100-1)/9 = 11）で分割して解を求める。各ステップのyの値を記録し、10x10の行列Aを作成（各列が各ステップのyベクトル）。
ベクトルbを[1,2,...,10]^Tとして、Ax=bをガウスの消去法で解く。

この流れで進めます。

まず、ルンゲクッタ法の実装。4次のルンゲクッタ法の場合、次のようにステップを進めます。各ステップで、k1, k2, k3, k4を計算し、次のyを求めます。

y_{n+1} = y_n + (h/6)(k1 + 2k2 + 2k3 + k4)

ここで、yがベクトルの場合、各kもベクトルとして計算されます。f(x, y)はベクトルを返す関数です。例えば、f(x, y) = -y （各要素を-1倍）のような単純な関数を仮定します。この場合、解はy(x) = y0 * exp(-(x - x0)) となりますが、これはテストしやすいです。

次に、行列Aの構築。例えば、初期値x0=1、y0は10次元ベクトル（例えば、最初の要素が1で他が0）とします。各ステップでxを11ずつ増やし（h=11）、10ステップでx=1 + 10*11 = 111となってしまいますが、これは問題文のx=100までの範囲を超えてしまいます。ステップ数とステップ幅を調整する必要があります。

正しくは、x=1からx=100までの区間を9分割して、10点（x=1, 1 + h, 1+2h, ...,1+9h=100）となるようにhを計算します。h = (100 -1)/9 = 99/9 = 11。つまり、ステップ幅h=11。これで10ステップ進むとx=1 + 9*11=100となります。各ステップでxを11ずつ増やし、各ステップでのyの値を記録します。その結果、10個のyベクトル（各10次元）が得られ、これを列として並べた10x10行列Aを作成します。

次にガウスの消去法の実装。行列Aとベクトルbから、連立方程式を解きます。部分ピボット選択付きのガウスの消去法を実装する必要があります。行列を上三角行列に変換し、後退代入を行います。

では、具体的なコードの流れを考えます。まず、numpyを使えるとのことなので、行列やベクトルの操作にnumpy.arrayを使用しますが、ルンゲクッタやガウス消去法のアルゴリズムは自作します。

コードの大まかな流れ：

f(x, y)を定義。例：f(x, y) = -y
初期条件y0を設定。x0=1、y0は10次元ベクトル（例：[1,0,...,0]）
ルンゲクッタ法でx=1からx=100までh=11で解を求め、各ステップのyを保存。
保存したyを列とする行列Aを作成。
ベクトルb = [1,2,...,10]^Tを作成。
ガウスの消去法でAx=bを解き、xを求める。

では、各部分を実装します。

まず、ルンゲクッタ法の関数。引数として、f、x0、y0、h、stepsを取る。各ステップでyを計算し、結果をリストに保存する。

def runge_kutta(f, x0, y0, h, steps):

y = y0.copy()

x = x0

result = [y0]

for _ in range(steps-1):

k1 = f(x, y)

k2 = f(x + h/2, y + (h/2)*k1)

k3 = f(x + h/2, y + (h/2)*k2)

k4 = f(x + h, y + h*k3)

y = y + (h/6)(k1 + 2k2 + 2*k3 + k4)

x += h

result.append(y.copy())

return np.array(result).T # 各列がステップの結果になるように転置

ただし、ここでyがnumpy配列であることを前提としています。例えば、y0が10次元のnumpyベクトルなら、k1などもベクトルとして計算されます。関数fはxとyを受け取り、dy/dxを返す。例としてfは：

def f(x, y):

return -y

初期値y0は、例えば：

y0 = np.zeros(10)

y0[0] = 1.0

次に、xの範囲を1から100まで、ステップ数10（h=11）で計算するため、steps=10。

result = runge_kutta(f, x0=1, y0=y0, h=11, steps=10)

すると、resultは10x10の行列（各列が各ステップのy）になります。これが行列Aです。

次に、ベクトルbを作成：

b = np.arange(1, 11)

そして、ガウスの消去法を実装します。ガウスの消去法の手順は、前進消去と後退代入です。部分ピボット選択を実装する必要があります。

def gauss_elimination(A, b):

n = len(b)

拡大係数行列を作成

aug = np.hstack([A, b.reshape(-1, 1)])

前進消去

for i in range(n):

部分ピボット選択

pivot_row = np.argmax(np.abs(aug[i:, i])) + i

aug[[i, pivot_row]] = aug[[pivot_row, i]]

ピボットが0なら特異行列

if aug[i, i] == 0:

raise ValueError("行列が特異です。")

ピボット行でi列目を1にするために、他の行を掃き出す

for j in range(i+1, n):

factor = aug[j, i] / aug[i, i]

aug[j, i:] -= factor * aug[i, i:]

後退代入

x = np.zeros(n)

for i in range(n-1, -1, -1):

x[i] = (aug[i, -1] - np.dot(aug[i, i+1:n], x[i+1:])) / aug[i, i]

return x

ただし、この実装では、拡大係数行列を作成し、前進消去で上三角行列に変換し、後退代入で解を求めています。部分ピボット選択を行い、行の入れ替えを行っています。

しかし、Aが10x10の行列で、bが10次元ベクトルなので、この関数にAとbを渡せば解xが得られます。

全体のコードをまとめると：

import numpy as np

微分方程式の右辺

def f(x, y):

return -y # 例として、dy/dx = -y

ルンゲクッタ法の実装

def runge_kutta(f, x0, y0, h, steps):

y = y0.copy()

x = x0

result = [y0.copy()]

for _ in range(steps - 1):

k1 = f(x, y)

k2 = f(x + h/2, y + (h/2)*k1)

k3 = f(x + h/2, y + (h/2)*k2)

k4 = f(x + h, y + h*k3)

y += (h/6) * (k1 + 2k2 + 2k3 + k4)

x += h

result.append(y.copy())

結果を転置して各列がステップのyになるように

return np.array(result).T

ガウスの消去法

def gauss_elimination(A, b):

n = len(b)

aug = np.hstack([A, b.reshape(-1, 1)])

前進消去

for i in range(n):

ピボット選択

pivot = np.argmax(np.abs(aug[i:, i])) + i

aug[[i, pivot]] = aug[[pivot, i]]

if aug[i, i] == 0:

raise ValueError("行列が特異です。")

掃き出し

for j in range(i+1, n):

factor = aug[j, i] / aug[i, i]

aug[j, i:] -= factor * aug[i, i:]

後退代入

x = np.zeros(n)

for i in range(n-1, -1, -1):

x[i] = (aug[i, -1] - aug[i, i+1:n] @ x[i+1:]) / aug[i, i]

return x

パラメータ設定

x0 = 1.0

y0 = np.zeros(10)

y0[0] = 1.0 # 初期値、例として最初の要素のみ1

h = (100 - 1) / 9 # 10ステップでx=100まで

steps = 10

ルンゲクッタで解を求め、行列Aを得る

A = runge_kutta(f, x0, y0, h, steps)

ベクトルbを作成

b = np.arange(1, 11)

ガウスの消去法でAx = bを解く

x = gauss_elimination(A, b)

print("Solution x:", x)

ただし、この例では微分方程式がdy/dx = -yであり、初期値y0が[1,0,...,0]なので、解はy(x) = y0 * exp(-(x - x0))となります。例えば、x=1でのyは[1,0,...,0]、x=12（次のステップ）でのyは各要素がexp(-11)倍されます。しかし、行列Aの各要素が非常に小さくなり、行列が特異に近くなる可能性があります。そのため、この例ではAx=bの解が存在しないか、数値的に不安定になる可能性があります。その場合、別のf(x, y)を定義する必要があります。

例えば、f(x, y) = y とすると、解は指数関数的に増加しますが、行列Aの要素が大きくなりすぎるかもしれません。または、線形代数の問題を解きやすくするため、Aが単位行列や対角行列になるような微分方程式を選ぶことも考えられます。例えば、dy/dx = 0とすると、yは定数になります。初期値y0が[1,0,...,0]なら、すべてのステップでyは同じになるので、行列Aはすべての列がy0となり、ランクが1の行列になります。この場合、Ax=bの解は存在しないか、不定になります。

そのため、微分方程式の右辺f(x, y)を調整して、行列Aが正則になるようにする必要があります。例えば、各要素が独立に変化するようなシステムを考える。例えば、dy_i/dx = -i * y_i のような場合、解はy_i(x) = y_i(x0) * exp(-i*(x - x0)) となります。異なる減衰率を持つため、行列Aは対角行列に近くなり、正則になる可能性が高いです。

そこで、f(x, y)を以下のように変更します：

def f(x, y):

n = len(y)

return -np.arange(1, n+1) * y # dy_i/dx = -i * y_i

初期値y0は、例えば全て1とします：

y0 = np.ones(10)

これにより、各成分y_iは異なる速度で減衰するため、行列Aは対角優位になり、正則になる可能性が高まります。

修正後のコード：

import numpy as np

def f(x, y):

return -np.arange(1, len(y)+1) * y # dy_i/dx = -i*y_i

def runge_kutta(f, x0, y0, h, steps):

y = y0.copy()

x = x0

result = [y0.copy()]

for _ in range(steps - 1):

k1 = f(x, y)

k2 = f(x + h/2, y + (h/2)*k1)

k3 = f(x + h/2, y + (h/2)*k2)

k4 = f(x + h, y + h*k3)

y += (h/6) * (k1 + 2k2 + 2k3 + k4)

x += h

result.append(y.copy())

return np.array(result).T

def gauss_elimination(A, b):

n = len(b)

aug = np.hstack([A, b.reshape(-1, 1)])

for i in range(n):

pivot = np.argmax(np.abs(aug[i:, i])) + i

aug[[i, pivot]] = aug[[pivot, i]]

if aug[i, i] == 0:

raise ValueError("Matrix is singular.")

for j in range(i+1, n):

factor = aug[j, i] / aug[i, i]

aug[j, i:] -= factor * aug[i, i:]

x = np.zeros(n)

for i in range(n-1, -1, -1):

x[i] = (aug[i, -1] - aug[i, i+1:n] @ x[i+1:]) / aug[i, i]

return x

x0 = 1.0

y0 = np.ones(10)

h = (100 - 1) / 9

steps = 10

A = runge_kutta(f, x0, y0, h, steps)

b = np.arange(1, 11)

try:

x = gauss_elimination(A, b)

print("Solution x:", x)

except ValueError as e:

print(e)

この場合、各y_iの減衰率が異なるため、行列Aは対角行列に近くなり、正則である可能性が高まります。これにより、ガウスの消去法で解が求まります。

テストとして、行列Aが対角行列であれば、解はbを各対角要素で割るだけで求まります。この例では、Aの各列は各ステップのyの値ですが、初期値が全て1で、dy_i/dx = -iy_iなので、x=1でのy_i(1) = 1。次のステップx=1 + h=1+11=12でのy_i(12) = exp(-i(12 -1 )) = exp(-i11). 同様に、各ステップでのy_i(x) = exp(-i(x -1 )). 行列Aの各列は、x=1, 12, 23, ...,100 でのy_iの値。したがって、行列Aの(j, k)要素は、y_j+1(x_k) = exp(- (j+1) * (x_k -1 )) となります。ここで、x_k = 1 + h*(k-1), h=11, k=1,...,10。例えば、k=1ならx=1, k=2ならx=12、など。そのため、行列Aの各列は異なる時間点でのyの値を表し、これが正則かどうかは指数関数の独立性によります。数学的には、異なる指数関数の組み合わせは線形独立であるため、行列Aは正則になる可能性が高いですが、数値的には非常に小さい値になるため、行列が悪条件になる可能性があります。

このコードを実行すると、行列Aが非常に小さい値（指数関数的に減少）を含むため、数値的に不安定になるかもしれません。その場合は、f(x, y)を調整するか、異なる初期条件を選ぶ必要があります。例えば、f(x, y) = y とすると、解は指数関数的に増加し、行列Aの要素が大きくなり、数値的により安定する可能性があります。

f(x, y)を以下のように変更：

def f(x, y):

return np.arange(1, len(y)+1) * y # dy_i/dx = i*y_i

この場合、解はy_i(x) = y_i(1) * exp(i*(x -1 )). 初期値y0が全て1なら、行列Aの各要素は大きな値になります。これにより、行列の条件数が改善される可能性があります。

この修正後、コードを実行すると、行列Aが大きくなるため、ガウスの消去法でも数値的に安定する可能性が高まります。しかし、指数関数的に増加するため、オーバーフローする可能性もあります。その場合は、hを小さくするか、ステップ数を減らす必要がありますが、問題文の指定によりx=1からx=100までを10ステップで進むためh=11である必要があります。

数値計算の観点から、このような大きなhを使うと、ルンゲクッタ法の誤差が大きくなる可能性がありますが、問題文では解法の正確さよりも、アルゴリズムの実装を求めていると考えられるので、この点は無視します。

以上を踏まえ、最終的なコードをまとめます。

import numpy as np

# 微分方程式の右辺を定義
def f(x, y):
    return np.arange(1, len(y)+1) * y  # dy_i/dx = i*y_i

# 4次のルンゲクッタ法を実装
def runge_kutta(f, x0, y0, h, steps):
    y = y0.copy()
    x = x0
    result = [y0.copy()]
    for _ in range(steps - 1):
        k1 = f(x, y)
        k2 = f(x + h/2, y + (h/2)*k1)
        k3 = f(x + h/2, y + (h/2)*k2)
        k4 = f(x + h, y + h*k3)
        y += (h/6) * (k1 + 2*k2 + 2*k3 + k4)
        x += h
        result.append(y.copy())
    return np.array(result).T  # 転置して各列がステップの結果に

# 部分ピボット選択付きガウスの消去法
def gauss_elimination(A, b):
    n = len(b)
    aug = np.hstack([A, b.reshape(-1, 1)]).astype(np.float64)
    
    for i in range(n):
        # ピボット選択
        pivot = np.argmax(np.abs(aug[i:, i])) + i
        aug[[i, pivot]] = aug[[pivot, i]]
        
        if aug[i, i] == 0:
            raise ValueError("Matrix is singular.")
        
        # 前進消去
        for j in range(i+1, n):
            factor = aug[j, i] / aug[i, i]
            aug[j, i:] -= factor * aug[i, i:]
    
    # 後退代入
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (aug[i, -1] - aug[i, i+1:n] @ x[i+1:]) / aug[i, i]
    return x

# パラメータ設定
x0 = 1.0
y0 = np.ones(10)  # 初期値（全て1）
h = (100 - x0) / 9  # ステップ幅（10ステップで100まで）
steps = 10

# 微分方程式を解いて行列Aを生成
A = runge_kutta(f, x0, y0, h, steps)

# ベクトルbを生成
b = np.arange(1, 11)

# 連立方程式を解く
try:
    x = gauss_elimination(A, b)
    print("解 x =", x)
except ValueError as e:
    print("エラー:", e)

このコードを走らせてみたところ、実行できてしまいました！o1でも苦戦した問題をあっさりコピペ一発で解きました。もう来年は「LLMには解けないけど人間には解ける」という問題は僕でも考え付けないかもですね。

何より、僕が張っていた罠である「何も考えないで適当に関数を定義して微分方程式を積分して係数行列Aに代入すると非正則になって解けない」という部分に言われなくても自分で気づいて関数を訂正し、しかも、その訂正はo1が考えたものよりスマートというところがですね、なんとも。