atividade3333

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define M_PI           3.14159265358979323846

double aceleracao(int tipo, double x, double k1x, double alfa, double f_externa, double gama, double v, double k1v)
{
    double a = 0;

    if(tipo == 1)
    {
        a = -alfa*(x + k1x/2) - gama*(v + k1v/2) + f_externa;
    }
    else if( tipo == 2)
    {
        a = -alfa*sin(x + k1x/2) - gama*(v + k1v/2)+ f_externa;
    }
    return a;
}

void euler(int tipo, double v, double x, double m, double l, double w0)
{
    double alfa = 0;
    double fx = 0;
    double gx = 0;
    double a = 0;
    double dt = 0.005;
    double epsilon = 0;
    double energia_mecanica = 0;
    double g = 9.8;
    double tempo = 0;
    double x_analitico = 0;
    double v_analitico = 0;
    double energia_analitica = 0;
    int i;
    char nome[30];
    char nome2[30];

    sprintf(nome, "euler-w%d.txt", tipo);
    sprintf(nome2, "diferenca-euler-w%d.txt", tipo);

    FILE *dados;
    FILE *diferenca;

    dados = fopen(nome, "w+");
    diferenca = fopen(nome2, "w+");

    for(i = 0; i <= 20000; i++)
    {
        tempo = i*dt;
        if(tipo == 1)
        {
            alfa = w0*w0;
            fx = x;
            gx = x*x/2;
            epsilon = v*v/2 + alfa*gx;
            energia_mecanica = epsilon*m;
        }
        else if(tipo == 2)
        {
            alfa = g/l;
            fx = sin(x);
            gx = 1.0 - cos(x);
            epsilon = v*v/2.0 + alfa*gx;
            energia_mecanica = epsilon*m*l*l;
        }
        fprintf(dados, " %f\t%f\t%f\t%f\n", x, v, energia_mecanica, tempo);
        fprintf(diferenca," %f\t%f\t%f\t%f\n", fabs(x - x_analitico), fabs(v - v_analitico), fabs(energia_mecanica - energia_analitica), tempo);
        a = -alfa*fx;
        x = x + v*dt;
        v = v + a*dt;

        x_analitico = 5*sin(w0*tempo)/w0;
        v_analitico = 5*cos(w0*tempo);

    }
    fclose(dados);
    fclose(dados);
}

void euler_cromer(int tipo, double v,double x, double m, double l, double w0)
{
    double alfa = 0;
    double fx = 0;
    double gx = 0;
    double a = 0;
    double dt = 0.001;
    double epsilon = 0;
    double energia_mecanica = 0;
    double g = 9.8;
    double tempo = 0;
    double gama = 0.0;
    double f_externa = 0;
    double wd = M_PI/3;
    int i;
    char nome[30];


    sprintf(nome, "euler_cromer-w%d.txt", tipo);

    FILE *dados;
    dados = fopen(nome, "w+");

    for(i = 0; i <= 20000; i++)
    {
        tempo = i*dt;
        if(tipo == 1)
        {
            alfa = w0*w0;
            fx = x;
            gx = x*x/2;
            epsilon = v*v/2 + alfa*gx;
            energia_mecanica = epsilon*m;
        }
        else if(tipo == 2)
        {
            alfa = g/l;
            fx = sin(x);
            gx = 1.0 - cos(x);
            epsilon = v*v/2.0 + alfa*gx;
            energia_mecanica = epsilon*m*l*l;
        }
        fprintf(dados, " %f\t%f\t%f\t%f\n", x, v, energia_mecanica, tempo);
        f_externa = 0; //0.05*cos(wd*tempo);
        a = -alfa*fx - gama*v + f_externa;
        v = v + a*dt;
        x = x + v*dt;

    }
    fclose(dados);
}

void kutta(int tipo, double v, double x, double m, double l, double w0)
{
    double alfa = 0;
    double fx = 0;
    double gx = 0;
    double a = 0;
    double dt = 0.001;
    double tempo = 0;
    double epsilon = 0;
    double energia_mecanica = 0;
    double g = 9.8;
    double k1v = 0;
    double k2v = 0;
    double k1x = 0;
    double k2x = 0;
    double x_antes = 0;
    double x_antes2 = 0;
    double gama = 0.2;
    double f_externa = 0;
    double wd = M_PI/10;
    int i;
    char nome[30];
    char nome2[30];

    sprintf(nome, "123-w%d.txt", tipo);
    sprintf(nome2, "amplitude-tempo.txt", tipo);

    FILE *dados;
    FILE *amplitude;
    amplitude = fopen(nome2, "w+");
    dados = fopen(nome, "w+");

    for(i = 0; i <= 20000; i++)
    {
        tempo = i*dt;
        if(tipo == 1)
        {
            alfa = w0*w0;
            fx = x;
            gx = x*x/2;
            epsilon = v*v/2 + alfa*gx;
            energia_mecanica = epsilon*m;
        }
        else if(tipo == 2)
        {
            alfa = g/l;
            fx = sin(x);
            gx = 1.0 - cos(x);
            epsilon = v*v/2.0 + alfa*gx;
            energia_mecanica = epsilon*m*l*l;
        }
        fprintf(dados, " %f\t%f\t%f\t%f\n", x, v, energia_mecanica, tempo);
        f_externa = 0.05*cos(wd*tempo);
        a = -alfa*fx - gama*v + f_externa;
        k1v = a*dt;
        k1x = v*dt;
        k2v = aceleracao(tipo, x, k1x, alfa, f_externa, gama, v, k1v)*dt;
        k2x =  (v + k1v/2)*dt;

        v = v + k2v;
        x_antes2 = x_antes;
        x_antes = x;
        x = x + k2x;
        if((x_antes - x_antes2 > 0 && x - x_antes < 0) || (x_antes - x_antes2 < 0 && x - x_antes > 0))
        {
            fprintf(amplitude, "%f\t %f\n", fabs(x_antes), tempo);
        }
    }
    fclose(dados);
    fclose(amplitude);
}

int main()
{
    double v = 0.2;
    double x = 0;
    double m = 1;
    double l = 1;
    double w0 = M_PI/2;

    euler(1, v, x, m, l, w0); // oscilador harmônico
    euler(2, v, x, m, l, w0); // pendulo simples
    euler_cromer(1, v, x, m, l, w0);
    euler_cromer(2, v, x, m, l, w0);
    kutta(1, v, x, m, l, w0);
    kutta(2, v, x, m, l, w0);

    return 0;
}