Metodo di Heun per risolvere equazioni differenziali del II ordine, sempre relativo ad un moto oscillatorio:



//Equazioni differenziali del II ordine
//Metodo di Heun
//Moto oscillatorio

#include < iostream >
#include < fstream >
#include < cmath >
#include "funz.h"
#include "TGraph.h"
#include "TApplication.h"
#include "TCanvas.h"

using namespace std;

//Condizioni iniziali

int main(int argc, char*argv[]) {

TApplication App("App", &argc, argv);
  double x[2], h, tmax;
  cerr << "Inserisci la posizione al tempo 0: ";
  cin  >> x[0];
  cerr << "Inserisci la velocita' al tempo 0: ";
  cin  >> x[1];
  cerr << "Integro fino al tempo: ";
  cin  >> tmax;
  cerr << "con passo:  ";
  cin  >> h;

//Calcolo della posizione al tempo finale

  ofstream f("oscill3.dat");
  double t=0;
  while (t<=tmax) {
    f << t << "\t" << x[0] << "\t" << x[1] << endl;
    heun(t,h,2,x,oscillatore);
  }
  f.close();

//Grafico tempo - posizione

  ifstream in1("oscill3.dat");
  double v,s;
  TGraph grafico1;
  int punto1=0;
  while (1) {
    in1 >> t >> s >> v;
    if ( in1.eof() ) break;
    grafico1.SetPoint(punto1,t,s);
    punto1++;
  }
 
  TCanvas *c1 = new TCanvas("c1", "Grafico posiz - tempo", 600, 400);
  c1->cd();
  grafico1.Draw("AL");
  in1.close();
 
//Grafico tempo - posizione + errori
 
  double fre=1. ; // frequenza
  double omega=2.*TMath::Pi()*fre;
  double x0=0., v0=1.; // condizioni iniziali
  double E = 0.5*(v0*v0+omega*omega*x0*x0);   // energia dell'oscillatore
  double A = sqrt(2.*E)/omega;   // ampiezza dell'oscillazione
  double Phi = atan2(v0,x0*omega);     // fase dell'oscillazione al tempo 0
 
  //double x,v;
  TGraph grafico2;
  int punto2=0;
 
  ifstream in2("oscill3.dat");
  while (1) {
    in2 >> t >> s >> v;
    if ( in2.eof() ) break;
    double err = fabs(s-A*cos(omega*t-Phi));
    grafico2.SetPoint(punto2,t,err);
    punto2++;
  }
  in2.close();

  TCanvas *c2 = new TCanvas("c2", "Grafico errori", 600, 400);
    c2->cd();
    grafico2.Draw("AL");

App.Run();
return 0;

}