pythonでdx/dt=-kxの式でフィッティングしたいのですがどうやってやればいいですか。
一次式の場合は、下記のコードで実行できたのですが常微分方程式のdx/dt=-kxでフィッティングしたいです。
よろしくお願いします。
x = np.array([20, 30, 40, 50, 60])
y = np.array([30, 38, 41, 49, 62])
p = np.polyfit(x, y, 1)
f = np.poly1d(p)
plt.scatter(x, y)
plt.plot(x, f(x))
せめて変数名は合わせてください。以下、$\frac{dy}{dx} = -ky$という問題だとして記載します。
まずは常微分方程式を解いてみましょう。$\frac{dy}{dx} = -ky$は1階線形常微分方程式の同次形といい、微分方程式の中では非常に簡単なタイプです。これを解いて$y =f(x)$の形にします。例えばこちらの解説を読んでみてください。
https://physnotes.jp/diffeq/1st-lde-gsol/
https://math0.pm.tokushima-u.ac.jp/~ohyama/lecture/diff_eq1/DE03.pdf (最初の3ページ)
ところで、提示されているxとyのサンプルは関数の形状に全然合っていないです。$y = f(x)$がどんな曲線になっていそうかイメージしてみてください。仮にdx=1,k=1と入れてみると、$dy = -y$となって、常にy方向の変化量がマイナスになります。つまりxが増えるにつれてyが減る単調減少関数であるはずです。
微分方程式の一般解$y = f(x)$が求められたら、次はフィッティングして未知の変数(kと積分定数C)を求めるわけですが、scipyのcurve_fit関数を使うのが良いかと思います。
カーブフィッティング手法 scipy.optimize.curve_fit の使い方を理解する
答えは載せますが、内容をよく理解するようにしてください。よろしくお願いします。
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
sample_x = np.array([0, 1, 2, 3, 4])
sample_y = np.array([1.0, 0.5, 0.13, 0.15, 0.04])
plt.scatter(sample_x, sample_y)
def f(x, C, k): # y = f(x)
return C * np.exp(-k * x)
# Fitting
C, k = curve_fit(f, sample_x, sample_y)[0]
print(f"y = {C:.6}e^(-{k:.6}x)")
# Plot
x = np.arange(0, 10, 0.01)
y = f(x, C, k)
plt.plot(x, y)
plt.show()
@xshangtiao
Questioner