method_Euler.cpp

#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

double B = 0.01, l = 1, E = 5, alpha = M_PI/4, m = 0.05, lambda = 0.1, g = 9.81, x0 = 0.1, dx0 = 0, t0 = 0, t1 = 20, eps = 0.1;

double func(double x, double dx)
{
    return (B*l*cos(alpha)/(m*lambda*x)) * (E + B*l*dx) - g * sin(alpha);
}

void Euler()
{
    ofstream fout;
    fout.open("Euler.txt");
    double dt = 1000000.;
    double x1 = x0, dx1 = dx0;
    double time = 0;
    int cnt = 1;
    while (time < (t1 - t0))
    {
        double x2, dx2;

        dx2 = dx1 + dt * func(x1, dx1);
        x2 = x1 + dt * dx2;

        while (abs(x2 - x1) > eps)
        {
            dx2 = dx1 + dt * func(x1, dx1);
            x2 = x1 + dt * dx2;
            dt /= 2;
        }

        fout << dx1 << " " << x1 << " " << time << endl;
        time += dt;
        dx1 = dx2;
        x1 = x2;
        cnt++;
    }
    fout << dx1 << " " << x1 << " " << time << endl;
    cout << time/cnt << endl;
    fout.close();
}

double dfunc(double x1, double dx1, double dt)
{
    double dx2;
    double k1, k2, k3, k4;

    k1 = func(x1, dx1);
    k2 = func(x1, dx1 + k1 * dt/2);
    k3 = func(x1, dx1 + k2 * dt/2);
    k4 = func(x1, dx1 + k3 * dt);
    dx2 = dx1 + (k1 + 2 * k2 + 2 * k3 + k4) * dt/6;

    return dx2;
}

void Runge()
{
    ofstream fout;
    fout.open("Runge.txt");
    double dt = 1000000.;
    double x1 = x0, dx1 = dx0;
    double time = 0;
    int cnt = 1;
    while (time < (t1 - t0))
    {
        double x2, dx2;
        double k1, k2, k3, k4;

        dx2 = dfunc(x1, dx1, dt);

        k1 = dfunc(x1, dx1, dt);
        k2 = dfunc(x1 + k1 * dt/2, dx1, dt);
        k3 = dfunc(x1 + k2 * dt/2, dx1, dt);
        k4 = dfunc(x1 + k3 * dt, dx1, dt);
        x2 = x1 + (k1 + 2 * k2 + 2 * k3 + k4) * dt/6;

        while (abs(x2 - x1) > eps)
        {
            dx2 = dfunc(x1, dx1, dt);

            k1 = dfunc(x1, dx1, dt);
            k2 = dfunc(x1 + k1 * dt/2, dx1, dt);
            k3 = dfunc(x1 + k2 * dt/2, dx1, dt);
            k4 = dfunc(x1 + k3 * dt, dx1, dt);
            x2 = x1 + (k1 + 2 * k2 + 2 * k3 + k4) * dt/6;

            dt /= 2;
        }

        fout << dx1 << " " << x1 << " " << time << endl;
        time += dt;
        dx1 = dx2;
        x1 = x2;
        cnt++;
    }
    fout << dx1 << " " << x1 << " " << time << endl;
    cout << time/cnt << endl;
    fout.close();
}

int main()
{
    int a;
    cin << a;
    if (a == 1)
    {
        Euler();
    }
    else
    {
        Runge();
    }
    return 0;
}

plot.py

import matplotlib.pyplot as plt

def read(file):
    time, f, df = [], [], []
    for line in file :
        line = line.split()
        time.append(float(line[2]))
        f.append(float(line[1]))
        df.append(float(line[0]))  
    return time, f, df

file1 = open("Runge.txt", "r")
file2 = open("Euler.txt", "r")
time_rk, f_rk, df_rk = read(file1)
time_eu, f_eu, df_eu = read(file2)

graph = plt.figure()
ax = graph.add_subplot(111)
ax.plot(time_eu, df_eu, "black")

graph = plt.figure()
ax = graph.add_subplot(111)
ax.plot(time_eu, f_eu, "black")

graph = plt.figure()
ax = graph.add_subplot(111)
ax.plot(time_rk, df_rk, "red")

graph = plt.figure()
ax = graph.add_subplot(111)
ax.plot(time_rk, f_rk, "red")

plt.show()

file1.close()
file2.close()