SHARE
TWEET

atividade3333

Tertius Dec 9th, 2019 90 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #define M_PI           3.14159265358979323846
  5.  
  6. double aceleracao(int tipo, double x, double k1x, double alfa, double f_externa, double gama, double v, double k1v)
  7. {
  8.     double a = 0;
  9.  
  10.     if(tipo == 1)
  11.     {
  12.         a = -alfa*(x + k1x/2) - gama*(v + k1v/2) + f_externa;
  13.     }
  14.     else if( tipo == 2)
  15.     {
  16.         a = -alfa*sin(x + k1x/2) - gama*(v + k1v/2)+ f_externa;
  17.     }
  18.     return a;
  19. }
  20.  
  21. void euler(int tipo, double v, double x, double m, double l, double w0)
  22. {
  23.     double alfa = 0;
  24.     double fx = 0;
  25.     double gx = 0;
  26.     double a = 0;
  27.     double dt = 0.005;
  28.     double epsilon = 0;
  29.     double energia_mecanica = 0;
  30.     double g = 9.8;
  31.     double tempo = 0;
  32.     double x_analitico = 0;
  33.     double v_analitico = 0;
  34.     double energia_analitica = 0;
  35.     int i;
  36.     char nome[30];
  37.     char nome2[30];
  38.  
  39.     sprintf(nome, "euler-w%d.txt", tipo);
  40.     sprintf(nome2, "diferenca-euler-w%d.txt", tipo);
  41.  
  42.     FILE *dados;
  43.     FILE *diferenca;
  44.  
  45.     dados = fopen(nome, "w+");
  46.     diferenca = fopen(nome2, "w+");
  47.  
  48.     for(i = 0; i <= 20000; i++)
  49.     {
  50.         tempo = i*dt;
  51.         if(tipo == 1)
  52.         {
  53.             alfa = w0*w0;
  54.             fx = x;
  55.             gx = x*x/2;
  56.             epsilon = v*v/2 + alfa*gx;
  57.             energia_mecanica = epsilon*m;
  58.         }
  59.         else if(tipo == 2)
  60.         {
  61.             alfa = g/l;
  62.             fx = sin(x);
  63.             gx = 1.0 - cos(x);
  64.             epsilon = v*v/2.0 + alfa*gx;
  65.             energia_mecanica = epsilon*m*l*l;
  66.         }
  67.         fprintf(dados, " %f\t%f\t%f\t%f\n", x, v, energia_mecanica, tempo);
  68.         fprintf(diferenca," %f\t%f\t%f\t%f\n", fabs(x - x_analitico), fabs(v - v_analitico), fabs(energia_mecanica - energia_analitica), tempo);
  69.         a = -alfa*fx;
  70.         x = x + v*dt;
  71.         v = v + a*dt;
  72.  
  73.         x_analitico = 5*sin(w0*tempo)/w0;
  74.         v_analitico = 5*cos(w0*tempo);
  75.  
  76.     }
  77.     fclose(dados);
  78.     fclose(dados);
  79. }
  80.  
  81. void euler_cromer(int tipo, double v,double x, double m, double l, double w0)
  82. {
  83.     double alfa = 0;
  84.     double fx = 0;
  85.     double gx = 0;
  86.     double a = 0;
  87.     double dt = 0.001;
  88.     double epsilon = 0;
  89.     double energia_mecanica = 0;
  90.     double g = 9.8;
  91.     double tempo = 0;
  92.     double gama = 0.0;
  93.     double f_externa = 0;
  94.     double wd = M_PI/3;
  95.     int i;
  96.     char nome[30];
  97.  
  98.  
  99.     sprintf(nome, "euler_cromer-w%d.txt", tipo);
  100.  
  101.     FILE *dados;
  102.     dados = fopen(nome, "w+");
  103.  
  104.     for(i = 0; i <= 20000; i++)
  105.     {
  106.         tempo = i*dt;
  107.         if(tipo == 1)
  108.         {
  109.             alfa = w0*w0;
  110.             fx = x;
  111.             gx = x*x/2;
  112.             epsilon = v*v/2 + alfa*gx;
  113.             energia_mecanica = epsilon*m;
  114.         }
  115.         else if(tipo == 2)
  116.         {
  117.             alfa = g/l;
  118.             fx = sin(x);
  119.             gx = 1.0 - cos(x);
  120.             epsilon = v*v/2.0 + alfa*gx;
  121.             energia_mecanica = epsilon*m*l*l;
  122.         }
  123.         fprintf(dados, " %f\t%f\t%f\t%f\n", x, v, energia_mecanica, tempo);
  124.         f_externa = 0; //0.05*cos(wd*tempo);
  125.         a = -alfa*fx - gama*v + f_externa;
  126.         v = v + a*dt;
  127.         x = x + v*dt;
  128.  
  129.     }
  130.     fclose(dados);
  131. }
  132.  
  133. void kutta(int tipo, double v, double x, double m, double l, double w0)
  134. {
  135.     double alfa = 0;
  136.     double fx = 0;
  137.     double gx = 0;
  138.     double a = 0;
  139.     double dt = 0.001;
  140.     double tempo = 0;
  141.     double epsilon = 0;
  142.     double energia_mecanica = 0;
  143.     double g = 9.8;
  144.     double k1v = 0;
  145.     double k2v = 0;
  146.     double k1x = 0;
  147.     double k2x = 0;
  148.     double x_antes = 0;
  149.     double x_antes2 = 0;
  150.     double gama = 0.2;
  151.     double f_externa = 0;
  152.     double wd = M_PI/10;
  153.     int i;
  154.     char nome[30];
  155.     char nome2[30];
  156.  
  157.     sprintf(nome, "123-w%d.txt", tipo);
  158.     sprintf(nome2, "amplitude-tempo.txt", tipo);
  159.  
  160.     FILE *dados;
  161.     FILE *amplitude;
  162.     amplitude = fopen(nome2, "w+");
  163.     dados = fopen(nome, "w+");
  164.  
  165.     for(i = 0; i <= 20000; i++)
  166.     {
  167.         tempo = i*dt;
  168.         if(tipo == 1)
  169.         {
  170.             alfa = w0*w0;
  171.             fx = x;
  172.             gx = x*x/2;
  173.             epsilon = v*v/2 + alfa*gx;
  174.             energia_mecanica = epsilon*m;
  175.         }
  176.         else if(tipo == 2)
  177.         {
  178.             alfa = g/l;
  179.             fx = sin(x);
  180.             gx = 1.0 - cos(x);
  181.             epsilon = v*v/2.0 + alfa*gx;
  182.             energia_mecanica = epsilon*m*l*l;
  183.         }
  184.         fprintf(dados, " %f\t%f\t%f\t%f\n", x, v, energia_mecanica, tempo);
  185.         f_externa = 0.05*cos(wd*tempo);
  186.         a = -alfa*fx - gama*v + f_externa;
  187.         k1v = a*dt;
  188.         k1x = v*dt;
  189.         k2v = aceleracao(tipo, x, k1x, alfa, f_externa, gama, v, k1v)*dt;
  190.         k2x =  (v + k1v/2)*dt;
  191.  
  192.         v = v + k2v;
  193.         x_antes2 = x_antes;
  194.         x_antes = x;
  195.         x = x + k2x;
  196.         if((x_antes - x_antes2 > 0 && x - x_antes < 0) || (x_antes - x_antes2 < 0 && x - x_antes > 0))
  197.         {
  198.             fprintf(amplitude, "%f\t %f\n", fabs(x_antes), tempo);
  199.         }
  200.     }
  201.     fclose(dados);
  202.     fclose(amplitude);
  203. }
  204.  
  205. int main()
  206. {
  207.     double v = 0.2;
  208.     double x = 0;
  209.     double m = 1;
  210.     double l = 1;
  211.     double w0 = M_PI/2;
  212.  
  213.     euler(1, v, x, m, l, w0); // oscilador harmônico
  214.     euler(2, v, x, m, l, w0); // pendulo simples
  215.     euler_cromer(1, v, x, m, l, w0);
  216.     euler_cromer(2, v, x, m, l, w0);
  217.     kutta(1, v, x, m, l, w0);
  218.     kutta(2, v, x, m, l, w0);
  219.  
  220.     return 0;
  221. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top