method_Euler_and_Runge-Kutty_two.cpp

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

using namespace std;

//parametres
double xo, yo, zo, b = 8/3, to, t1, eps;
int sigma = 10;

double F(double x, double y) // x' = F(x,y)
{
    return sigma*(y - x);
}

double G(double x, double y, double z, double r) // y' = G(x, y, z, r)
{
    return - x*z + r*x - y;
}

double H(double x, double y, double z) // z' = H(x, y, z)
{
    return x*y - b*z;
}

void Euler(double r) // Euler's method
{
    ofstream fout;
    fout.open("Euler.txt");
    double dt = 1000000.;
    double x1 = xo, y1 = yo, z1 = zo;
    double cnt_dt = 0;
    int i = 1;
    while (cnt_dt < (t1 - to))
    {
        double x2, y2, z2;
        x2 = x1 + dt * F(x1, y1);
        y2 = y1 + dt * G(x1, y1, z1, r);
        z2 = z1 + dt * H(x1, y1, z1);
        while (((abs(x1 - x2) > eps) && (abs(y1 - y2) > eps) && (abs(z1 - z2) > eps)) || (abs(dt - 1000000.) < 0.0001))
        {
            x2 = x1 + dt * F(x1, y1);
            y2 = y1 + dt * G(x1, y1, z1, r);
            z2 = z1 + dt * H(x1, y1, z1);
            dt /= 2;
        }
        fout << x1 << " " << y1 << " " << z1 << " " << cnt_dt << endl;
        cnt_dt += dt;
        x1 = x2;
        y1 = y2;
        z1 = z2;
        i++;
    }
    fout << x1 << " " << y1 << " " << z1 << " " << cnt_dt << endl;
    cout << cnt_dt/i;
    fout.close();
}

void Runge(double r) //Runge-Kutta method
{
    ofstream fout;
    fout.open("Runge.txt");
    double dt = 1000000.;
    double x1 = xo, y1 = yo, z1 = zo;
    double cnt_dt = 0;
    int i = 1;
    while (cnt_dt < (t1 - to))
    {
        double x2, y2, z2, k1, k2, k3, k4;
        k1 = F(x1, y1);
        k2 = F(x1+k1*dt/2, y1+k1*dt/2);
        k3 = F(x1+k2*dt/2, y1+k2*dt/2);
        k4 = F(x1+k3*dt, y1+k3*dt);
        x2 = x1 + ( k1 + 2*k2 + 2*k3 +k4 )*dt/6;

        k1 = G(x1, y1, z1, r);
        k2 = G(x1 + k1*dt/2, y1 + k1*dt/2, z1 + k1*dt/2, r);
        k3 = G(x1 + k2*dt/2, y1 + k2*dt/2, z1 + k2*dt/2, r);
        k4 = G(x1 + k3*dt, y1 + k3*dt, z1 + k3*dt, r);
        y2 = y1 + ( k1 + 2*k2 + 2*k3 + k4 )*dt/6;

        k1 = H(x1, y1, z1);
        k2 = H(x1 + k1*dt/2, y1 + k1*dt/2, z1 + k1*dt/2);
        k3 = H(x1 + k2*dt/2, y1 + k2*dt/2, z1 + k2*dt/2);
        k4 = H(x1 + k3*dt, y1 + k3*dt, z1 + k3*dt);
        z2 = z1 + ( k1 + 2*k2 + 2*k3 + k4 )*dt/6;

        while (((abs(x1 - x2) > eps) && (abs(y1 - y2) > eps) && (abs(z1 - z2) > eps)) || (abs(dt - 1000000.) < 0.0001))
        {
            k1 = F(x1, y1);
            k2 = F(x1+k1*dt/2, y1+k1*dt/2);
            k3 = F(x1+k2*dt/2, y1+k2*dt/2);
            k4 = F(x1+k3*dt, y1+k3*dt);
            x2 = x1 + ( k1 + 2*k2 + 2*k3 +k4 )*dt/6;

            k1 = G(x1, y1, z1, r);
            k2 = G(x1 + k1*dt/2, y1 + k1*dt/2, z1 + k1*dt/2, r);
            k3 = G(x1 + k2*dt/2, y1 + k2*dt/2, z1 + k2*dt/2, r);
            k4 = G(x1 + k3*dt, y1 + k3*dt, z1 + k3*dt, r);
            y2 = y1 + ( k1 + 2*k2 + 2*k3 + k4 )*dt/6;

            k1 = H(x1, y1, z1);
            k2 = H(x1 + k1*dt/2, y1 + k1*dt/2, z1 + k1*dt/2);
            k3 = H(x1 + k2*dt/2, y1 + k2*dt/2, z1 + k2*dt/2);
            k4 = H(x1 + k3*dt, y1 + k3*dt, z1 + k3*dt);
            z2 = z1 + ( k1 + 2*k2 + 2*k3 + k4 )*dt/6;
            
            dt /= 2;
        }
        fout << x1 << " " << y1 << " " << z1 << " " << cnt_dt << endl;
        cnt_dt += dt;
        x1 = x2;
        y1 = y2;
        z1 = z2;
        i++;
    }
    fout << x1 << " " << y1 << " " << z1 << " " << cnt_dt << endl;
    cout << cnt_dt/i;
    fout.close();
}

int main()
{
    double r;
    cin >> r;
    cin >> xo >> yo >> zo >> to >> t1 >> eps; //ввод начальных координат, промежутка времени, погрешности
    cout << "What's method? 1 is Euler. 2 is Runge." << endl;
    int kek = 0;
    cin >> kek;
    if (kek == 1) {Euler(r);}
    if (kek == 2) {Runge(r);}
    return 0;
}

plots.py

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d, Axes3D

def read(file):
    x, y, z, t = [], [], [], []
    for line in file :
        line = line.split()
        x.append(float(line[0]))
        y.append(float(line[1]))
        z.append(float(line[2]))  
        t.append(float(line[3]))
    return x, y, z, t

file = open("Runge.txt", "r")
file1 = open("Euler.txt", "r")
x, y, z, t = read(file)
x1, y1, z1, t1 = read(file1)

fig = plt.figure()
ax = fig.gca(projection="3d")
ax.plot(x, y, z, "red")
ax.plot(x1, y1, z1, "black")

graph = plt.figure()
ax = graph.add_subplot(111)
ax.plot(t1, x1, "black")
ax.plot(t, x, "red")

plt.show()