這是一系列的文章，它們的內容會從簡單到複雜。我打算首先介紹常微分方程組作圖的核心：“龍格-庫塔（Runge-Kutta)法”，然後從一維線性諧振子開始，到三體問題。我會把我寫的程式碼和影象貼出來供參考。
對於方程
，
，其中
，有
。
那麼
。
如果我們想提高精度y我們可以有一下形式
這個處理過程告訴我們，如果設法在
![\left[ x_{n},x_{n+1}\right]](https://img.bearask.com/upload/website_attach/EyJY1YWiwGxg2jgeOb4xOZ=fG_B+kaluB9xYHEkd8e4T3kdKqEG7NFNZso8p_bGNy70q00ta29sZKj3ygliMP1=M6KxJ8Y2+ooXLHEyQN9sy6XoJm.jpeg)
內多預報幾個點的斜率值，然後將它們加權平均作為平均斜率
，則有可能構造出更高精度的計算格式，這就是龍格-庫塔（Runge-Kutta）法。
我們一般用四階龍格-庫塔方法來處理常微分方程
四階龍格-庫塔法
經典的四階龍格-庫塔格式每一步需要計算四次函式值
，它具有四階精度。
例項
對於方程
，其中
這裡，我們可以輕易求得
的解析解，即
。而
。
我們用上面的龍格-庫塔公式驗證一下，取步長
輸入程式碼
#include
FILE
*
fp
=
fopen
（
“ex的值。dat”
，
“w”
）；
double
func
（
double
x
）
{
return
x
；
}
void
tworder
（
double
x_0
，
double
h
，
double
t_0
，
double
t_n
）
{
double
n
=
（
t_n
-
t_0
）
/
h
；
double
t
=
t_0
，
x
=
x_0
；
double
k1
，
k2
，
k3
，
k4
；
int
i
；
printf
（
“%f
\t
%f
\n
”
，
t
，
x
）；
fprintf
（
fp
，
“%f
\t
%f
\n
”
，
t
，
x
）；
for
（
i
=
1
；
i
<=
n
；
i
++
）
{
k1
=
func
（
x
）；
k2
=
func
（
x
+
k1
*
h
/
2
）；
k3
=
func
（
x
+
k2
*
h
/
2
）；
k4
=
func
（
x
+
k3
*
h
）；
x
=
x
+
（
k1
+
2
*
k2
+
2
*
k3
+
k4
）
*
h
/
6
；
t
=
t
+
h
；
printf
（
“%f
\t
%f
\n
”
，
t
，
x
）；
fprintf
（
fp
，
“%f
\t
%f
\n
”
，
t
，
x
）；
}
}
int
main
（）
{
double
h
=
0。01
；
double
x_0
=
1
；
double
t_0
=
0
，
t_n
=
10
；
tworder
（
x_0
，
h
，
t_0
，
t_n
）；
fclose
（
fp
）；
return
0
；
}
執行後得到
。
用龍格-庫塔法計算的結果並不是絕對精確。但是，只要時間
取值不要太大，誤差的範圍仍然是可以接受的。
