C++で常微分方程式(オイラー法)

をオイラー法で近似する

#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

// 元の関数
double f(double x);
// 導関数
double fd(double x);
// Euler法
void euler(double x, double y, int n, double max);

int main(void)
{
    double x = 0.0;
    double y = f(x);
    double max = 4;

    for (int n = 4; n <= 128; n *= 2) {
        cout << "分割数 = ";
        cout << n << endl;
        euler(x, y, n, max);
    }
    return 0;
}

// 元の関数
double f(double x)
{
    return x - pow(x,3) / (3 * 2) + pow(x,5) / (5 * 4 * 3 * 2);
}
// 導関数
double fd(double x)
{
    return 1 - pow(x,2) / 2 + pow(x,4) / (4 * 3 * 2);
}
// Euler法
void euler(double x, double y, int n, double max)
{
    double h = max / n;

    for (int i = 1; i <= n; i++) {
        y = y + h * fd(x);
        x = x + h;
        double z = f(x);

        // 元の関数と比較
        cout << setw(8) << fixed << setprecision(5) << x     << "\t";
        cout << setw(8) << fixed << setprecision(5) << y     << "\t";
        cout << setw(8) << fixed << setprecision(5) << z     << "\t";
        cout << setw(8) << fixed << setprecision(5) << y - z << endl;
    }
}

分割数 = 4
 1.00000	 1.00000	 0.84167	 0.15833
 2.00000	 1.54167	 0.93333	 0.60833
 3.00000	 1.20833	 0.52500	 0.68333
 4.00000	 1.08333	 1.86667	-0.78333
分割数 = 8
 0.50000	 0.50000	 0.47943	 0.02057
 1.00000	 0.93880	 0.84167	 0.09714
 1.50000	 1.20964	 1.00078	 0.20885
 2.00000	 1.25260	 0.93333	 0.31927
 2.50000	 1.08594	 0.70964	 0.37630
 3.00000	 0.83724	 0.52500	 0.31224
 3.50000	 0.77474	 0.73099	 0.04375
 4.00000	 1.33854	 1.86667	-0.52813
分割数 = 16
 0.25000	 0.25000	 0.24740	 0.00260
 0.50000	 0.49223	 0.47943	 0.01280
 0.75000	 0.71163	 0.68167	 0.02996
 1.00000	 0.89461	 0.84167	 0.05295
 1.25000	 1.03003	 0.94991	 0.08012
 1.50000	 1.11015	 1.00078	 0.10937
 1.75000	 1.13163	 0.99355	 0.13809
 2.00000	 1.09652	 0.93333	 0.16318
 2.25000	 1.01318	 0.83210	 0.18108
 2.50000	 0.89734	 0.70964	 0.18770
 2.75000	 0.77299	 0.59449	 0.17850
 3.00000	 0.67342	 0.52500	 0.14842
 3.25000	 0.64217	 0.55024	 0.09194
 3.50000	 0.73401	 0.73099	 0.00302
 3.75000	 1.01591	 1.14075	-0.12484
 4.00000	 1.56803	 1.86667	-0.29863
分割数 = 32
 0.12500	 0.12500	 0.12467	 0.00033
 0.25000	 0.24902	 0.24740	 0.00162
 0.37500	 0.37014	 0.36627	 0.00387
 0.50000	 0.48645	 0.47943	 0.00703
 0.62500	 0.59615	 0.58510	 0.01105
 0.75000	 0.69753	 0.68167	 0.01587
 0.87500	 0.78903	 0.76762	 0.02141
 1.00000	 0.86923	 0.84167	 0.02756
以下略

参考文献