false.cpp

#include <iostream>
#include <fstream>
#include <math.h>
#include <vector>

using namespace std;

double Uo, w, B, S, R, J, to, t1, fo, h;
/*
#define Uo 10
#define w 10
#define B 1
#define S 0.01
#define R 10
#define J 1
#define to 0
#define t1 3
#define fo 0.5
#define h 0.001
*/
#define n (int)(t1 - to)/h

double func(double f, double t,  double df)
{
    return (Uo*cos(w*t) + B*S*sin(f)*df)*S*B*sin(f)/(R*J);
}

double dfunc(double f, double t,  double df, double dt)
{
    double k1, k2, k3, k4;
    k1 = func(f, t, df);
    k2 = func(f + dt*k1/2., t + dt/2., df + dt*k1/2.);
    k3 = func(f + dt*k2/2., t + dt/2., df + dt*k2/2.);
    k4 = func(f + dt*k3, t + dt, df + dt*k3);
    return df + (k1 + 2*k2 + 2*k3 + k4)*dt/6.;
}

double g(double f, double t,  double df, double dt)
{
    double k1, k2, k3, k4;
    k1 = dfunc(f, t, df, dt);
    k2 = dfunc(f + dt*k1/2., t + dt/2., df + dt*k1/2., dt);
    k3 = dfunc(f + dt*k2/2., t + dt/2., df + dt*k2/2., dt);
    k4 = dfunc(f + dt*k3, t + dt, df + dt*k3, dt);
    return f + (k1 + 2*k2 + 2*k3 + k4)*dt/6.;
}

void Runge()
{
    ofstream fout;
    fout.open("RK.txt");
    vector<double> arr_df, arr_f, time;
    arr_df.push_back(0);
    arr_f.push_back(fo);
    time.push_back(to);
    double dt = h;
    for (int i = 1; i < n; i++)
    {
        arr_df.push_back( dfunc(arr_f[i-1], time[i-1], arr_df[i-1], dt) );
        arr_f.push_back( g(arr_f[i-1], time[i-1], arr_df[i-1], dt) );
        time.push_back(time[i-1] + dt);
        fout << time[i-1] << " " << arr_f[i-1] << " " << arr_df[i-1] << endl;
    }
    fout << time[n-1] << " " << arr_f[n-1] << " " << arr_df[n-1] << endl;
    fout.close();
}

void Euler()
{
    ofstream fout;
    fout.open("Euler.txt");
    vector<double> arr_df, arr_f, time;
    arr_df.push_back(0);
    arr_f.push_back(fo);
    time.push_back(to);
    double dt = h;
    for (int i = 1; i < n; i++)
    {
        arr_df.push_back( arr_df[i-1]  + dt*func(arr_f[i-1], time[i-1], arr_df[i-1]) );
        arr_f.push_back( arr_f[i-1] + dt*arr_df[i] );
        time.push_back(time[i-1] + dt);
        fout << time[i-1] << " " << arr_f[i-1] << " " << arr_df[i-1] << endl;
    }
    fout << time[n-1] << " " << arr_f[n-1] << " " << arr_df[n-1] << endl;
    fout.close();
}

int main()
{
    cin >> Uo >> w >> B >> S >> R >> J >> to >> t1 >> fo >> h;
    Runge();
    Euler();
    return 0;
}

method_Euler_and_Runge-Kutty.cpp

#include <iostream>
#include <fstream>
#include <math.h>
#include <vector>

using namespace std;

double Uo, w, B, S, R, J, to, t1, fo;
int n;

/*
#define Uo 10
#define w 10
#define B 1
#define S 0.01
#define R 10
#define J 1
#define to 0
#define t1 3
#define fo 0.5
#define n 100000
*/

double func(double f, double t,  double df)
{
    return (Uo*cos(w*t) + B*S*sin(f)*df)*S*B*sin(f)/(R*J);
}

double dfunc(double f, double t,  double df, double dt)
{
    double k1, k2, k3, k4;
    k1 = func(f, t, df);
    k2 = func(f + dt*k1/2., t + dt/2., df + dt*k1/2.);
    k3 = func(f + dt*k2/2., t + dt/2., df + dt*k2/2.);
    k4 = func(f + dt*k3, t + dt, df + dt*k3);
    return df + (k1 + 2*k2 + 2*k3 + k4)*dt/6.;
}

double g(double f, double t,  double df, double dt)
{
    double k1, k2, k3, k4;
    k1 = dfunc(f, t, df, dt);
    k2 = dfunc(f + dt*k1/2., t + dt/2., df + dt*k1/2., dt);
    k3 = dfunc(f + dt*k2/2., t + dt/2., df + dt*k2/2., dt);
    k4 = dfunc(f + dt*k3, t + dt, df + dt*k3, dt);
    return f + (k1 + 2*k2 + 2*k3 + k4)*dt/6.;
}

void Runge()
{
    ofstream fout;
    fout.open("RK.txt");
    vector<double> arr_df, arr_f, time;
    arr_df.push_back(0);
    arr_f.push_back(fo);
    time.push_back(to);
    double dt = double (t1 - to)/n;
    for (int i = 1; i < n; i++)
    {
        arr_df.push_back( dfunc(arr_f[i-1], time[i-1], arr_df[i-1], dt) );
        arr_f.push_back( g(arr_f[i-1], time[i-1], arr_df[i-1], dt) );
        time.push_back(time[i-1] + dt);
        fout << time[i-1] << " " << arr_f[i-1] << " " << arr_df[i-1] << endl;
    }
    fout << time[n-1] << " " << arr_f[n-1] << " " << arr_df[n-1] << endl;
    fout.close();
}

void Euler()
{
    ofstream fout;
    fout.open("Euler.txt");
    vector<double> arr_df, arr_f, time;
    arr_df.push_back(0);
    arr_f.push_back(fo);
    time.push_back(to);
    double dt = double (t1 - to)/n;
    for (int i = 1; i < n; i++)
    {
        arr_df.push_back( arr_df[i-1]  + dt*func(arr_f[i-1], time[i-1], arr_df[i-1]) );
        arr_f.push_back( arr_f[i-1] + dt*arr_df[i] );
        time.push_back(time[i-1] + dt);
        fout << time[i-1] << " " << arr_f[i-1] << " " << arr_df[i-1] << endl;
    }
    fout << time[n-1] << " " << arr_f[n-1] << " " << arr_df[n-1] << endl;
    fout.close();
}

int main()
{
    cin >> Uo >> w >> B >> S >> R >> J >> to >> t1 >> fo >> n;
    Runge();
    Euler();
    return 0;
}

plots.py

import matplotlib.pyplot as plt

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

file1 = open("RK.txt", "r")
file2 = open("Euler.txt", "r")
time, f_rk, df_rk = read(file1)
time, f_eu, df_eu = read(file2)
    
graph = plt.figure()
ax = graph.add_subplot(111)
ax.plot(time, df_rk, "red")
ax.plot(time, df_eu, "black")

graph = plt.figure()
ax = graph.add_subplot(111)
ax.plot(time, f_rk, "red")
ax.plot(time, f_eu, "black")

plt.show()

file1.close()
file2.close()