Euler_Runge.cpp

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

using namespace std;

double g = 9.8, R, p0, p, z0, v0 = 0, t, eps;

double func(double z)
{
    if (z >= R)
    {
        return g - g * (p0/p);
    }
    if (z <= -R)
    {
        return g;
    }
    return g - g * (p0/p) * (1/2 + 3*z/(4*R) - (1/4)*pow(z/R, 3));
}

void Euler()
{
    ofstream fout;
    fout.open("Euler.txt");
    double dt = 1000000.;
    double z1 = z0, v1 = v0;
    double time = 0;
    int cnt = 1;
    while (time < t)
    {
        double z2, v2;
        v2 = v1 + func(z1) * dt;
        z2 = z1 + v2 * dt;
        while (abs(z1 - z2) > eps)
        {
            dt /= 2;
            v2 = v1 + func(z1) * dt;
            z2 = z1 + v2 * dt;
        }
        fout << z1 << " " << v1 << " " << time << endl;
        time += dt;
        v1 = v2;
        z1 = z2;
        cnt++;
    }

    fout << z1 << " " << v1 << " " << time << endl;
    cout << time/cnt << endl;
    fout.close();
}

double dfunc(double z1, double v1, double dt)
{
    double v2;
    double k1, k2, k3, k4;

    k1 = func(z1);
    k2 = func(z1);
    k3 = func(z1);
    k4 = func(z1);
    v2 = v1 + (k1 + 2 * k2 + 2 * k3 + k4) * dt/6;

    return v2;
}

void Runge()
{
    ofstream fout;
    fout.open("Runge.txt");
    double dt = 1000.;
    double z1 = z0, v1 = v0;
    double time = 0;
    int cnt = 1;
    while (time < t)
    {
        double z2, v2;
        double k1, k2, k3, k4;

        v2 = dfunc(z1, v1, dt);
        k1 = dfunc(z1, v1, dt);
        k2 = dfunc(z1 + dt*k1/2, v1, dt);
        k3 = dfunc(z1 + dt*k2/2, v1, dt);
        k4 = dfunc(z1 + k3*dt, v1, dt);
        z2 = z1 + (k1 + 2 * k2 + 2 * k3 + k4)*dt/6;
        while (abs(z1 - z2) > eps)
        {
            dt /= 2;
            v2 = dfunc(z1, v1, dt);
            k1 = dfunc(z1, v1, dt);
            k2 = dfunc(z1 + dt*k1/2, v1, dt);
            k3 = dfunc(z1 + dt*k2/2, v1, dt);
            k4 = dfunc(z1 + k3*dt, v1, dt);
            z2 = z1 + (k1 + 2 * k2 + 2 * k3 + k4)*dt/6;
        }
        fout << z1 << " " << v1 << " " << time << endl;
        time += dt;
        v1 = v2;
        z1 = z2;
        cnt++;
    }

    fout << z1 << " " << v1 << " " << time << endl;
    cout << time/cnt << endl;
    fout.close();
}

int main()
{
    cout << "R p0 p z0 t eps" << endl;
    cin >> R >> p0 >> p >> z0 >> t >> eps;
    Euler();
    Runge();
    /*int a;
    cin >> a;
    if (a == 1) {Euler();}
    if (a == 2) {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[0]))
        df.append(float(line[1]))  
    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()