дифференциальные уравнения в с ++, задачи с методом Рунге Кутты

Я написал программу на C ++, которая должна решать дифференциальные уравнения. Проблема в том, что, похоже, он плохо работает с ROOT. Он компилируется нормально, но когда я выполняю, это то, что я получаю:

*** Break *** segmentation violation===========================================================
There was a crash.
This is the entire stack trace of all threads:
===========================================================
#0  0x00007fc28193984a in __GI___waitpid (pid=7730, stat_loc=stat_loc
entry=0x7ffffe4ae000, options=options
entry=0) at ../sysdeps/unix/sysv/linux/waitpid.c:31
#1  0x00007fc2818b2ffb in do_system (line=<optimized out>) at ../sysdeps/posix/system.c:148
#2  0x00007fc2831d0954 in TUnixSystem::StackTrace() () from /usr/lib/root/libCore.so
#3  0x00007fc2831d29ec in TUnixSystem::DispatchSignals(ESignals) () from /usr/lib/root/libCore.so
#4  <signal handler called>
#5  0x0000000000405a8a in Runge_Kutta::Passo(double, VettoreLineare&, double) ()
#6  0x0000000000403b8a in main ()
===========================================================The lines below might hint at the cause of the crash.
If they do not help you then please submit a bug report at
http://root.cern.ch/bugs. Please post the ENTIRE stack trace
from above as an attachment in addition to anything else
that might help us fixing this issue.
===========================================================
#5  0x0000000000405a8a in Runge_Kutta::Passo(double, VettoreLineare&, double) ()
#6  0x0000000000403b8a in main ()
===========================================================

Это моя программа

equazione_differenziale.c

#include "equazione_differenziale.h"
EqDifferenzialeBase :: EqDifferenzialeBase (FunzioneBase* f) {

_f=f;

};

Eulero :: Eulero (FunzioneBase*f) : EqDifferenzialeBase(f) {};

Runge_Kutta :: Runge_Kutta (FunzioneBase* f) : EqDifferenzialeBase(f) {};

VettoreLineare Eulero :: Passo (double t, VettoreLineare& x, double h) {

VettoreLineare vec(x.GetN());
vec=x+_f->Eval(t,x)*h;

return vec;

};

Protone :: Protone (double m, double q, double E0, double f, double lambda){
_m=m;
_q=q;
_E0=E0;
_f=f;
_lambda=lambda;
};VettoreLineare Protone::Eval(double t, const VettoreLineare& v) const{

VettoreLineare y(v.GetN());

for(int i=0; i<v.GetN()/2; i++){
y.SetComponent(i, v.GetComponent(v.GetN()/2+i));
y.SetComponent(i+v.GetN()/2, (-1.)*(_q/_m)*_E0*sin(2*M_PI*(v.GetComponent(i)/_lambda)-2*M_PI*_f*t));
};

return y;

};

VettoreLineare Runge_Kutta::Passo(double t, VettoreLineare& v, double h){

VettoreLineare k1=_f->Eval(t,v);
VettoreLineare k2=_f->Eval(t+h/2.,v+k1*(h/2.));
VettoreLineare k3=_f->Eval(t+h/2.,v+k2*(h/2.));
VettoreLineare k4=_f->Eval(t+h,v+k3*h);
VettoreLineare y=v+(k1+k2*2.+k3*2.+k4)*(h/6.);

return y;

};

equazione_differenziale.h

#ifndef equazione_differenziale_h_
#define equazione_differenziale_h_
#include "Vettore.h"#include <iostream>
#include <cmath>

class FunzioneBase {

public:
virtual VettoreLineare Eval(double t, const VettoreLineare& v) const=0;

};

class Protone: public FunzioneBase {

private:
double _m,_q,_E0,_f,_lambda;

public:
Protone(double m, double q, double E0, double f, double lambda);
virtual VettoreLineare Eval(double t, const VettoreLineare& v) const;};class EqDifferenzialeBase {

protected:
FunzioneBase* _f;

public:
EqDifferenzialeBase (FunzioneBase* f);
virtual VettoreLineare Passo (double t, VettoreLineare& x, double h)=0;
};

class Eulero : public EqDifferenzialeBase {

public:
Eulero (FunzioneBase* f);
virtual VettoreLineare Passo (double t, VettoreLineare& x, double h);

};

class Runge_Kutta: public EqDifferenzialeBase {

protected:
FunzioneBase* _f;

public:
Runge_Kutta (FunzioneBase* f);
virtual VettoreLineare Passo(double t, VettoreLineare& v, double h);
};#endif

Vettore.h

#ifndef vettore_h_
#define vettore_h_

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

class Vettore {

protected:
unsigned int _N;
double * _v;
void Quicksort(unsigned int primo, unsigned int ultimo);
void Scambia (int a, int b);

public:
Vettore ();
Vettore (int N);
Vettore (int N, char* nomefile);
Vettore (const Vettore& v);
virtual void SetComponent (int i, double x);
void AddComponent (double x);
double GetComponent (int i) const;
void Print () const;
void Print (char* nomefile) const;
void Sort();
virtual int GetN() const;
Vettore& operator=(const Vettore & vetty);
~Vettore();

};

class VettoreLineare : public Vettore {

protected:public:
VettoreLineare () : Vettore() {};
VettoreLineare (int N) : Vettore(N) {};
VettoreLineare (int N, char* nomefile) : Vettore(N, nomefile) {};
VettoreLineare (const Vettore& v) : Vettore(v) {};
VettoreLineare operator+(const VettoreLineare& v);
VettoreLineare operator*(double lambda);
VettoreLineare& operator=(const VettoreLineare& v);
virtual int GetN() const;
virtual void SetComponent(int i, double x);

};

Vettore.c

#include "Vettore.h"
//Default Constructor

Vettore :: Vettore () {

_N=0;
_v=NULL;

};

//N Constructor

Vettore :: Vettore (int N) {

_N=N;
_v=new double [_N];

for (int i=0; i<_N; i++)
_v[i]=0;

};

//N file-taken constructor

Vettore :: Vettore (int N, char* nomefile) {

_N=N;
_v=new double [_N];

std::ifstream input;
input.open(nomefile);

double dato;

input>>dato;

for(int i=0; i<N; i++){
_v[i]=dato;
input>>dato;
};

input.close();

};

//Copyconstructor

Vettore :: Vettore (const Vettore& V) {

_N=V.GetN();

_v=new double [_N];

for(int i=0; i<_N; i++)
_v[i]=V.GetComponent(i);};

//Destructor

Vettore::~Vettore(){

delete[] _v;

};

//Set Component

void Vettore :: SetComponent (int i, double x) {

if (i>_N) {
std::cout<<"errore!"<<std::endl;
return ;
};_v[i]=x;

};

//Get Component

double Vettore :: GetComponent (int i) const {

if (i>_N){
std::cout<<"errore!"<<std::endl;
return 0;
};

return _v[i];

};

//Add Component (aggiunge il valore desiderato nella coda del vettore)

void Vettore :: AddComponent (double x) {

double* a=new double [_N+1];

for(int i=0; i<_N; i++)
a[i]=_v[i];

a[_N]=x;
delete [] _v;
_v=a;

_N=_N+1;};

//Print

void Vettore :: Print () const {

std::cout<<"Il vettore ha: "<<_N<<" componenti."<<std::endl;

for(int i=0; i<_N; i++)
std::cout<<_v[i]<<std::endl;

};

//Stampa su file

void Vettore :: Print (char* nomefile) const {
std::ofstream output;
output.open(nomefile);

output<<_N;

for(int i=0; i<_N; i++)
output<<_v[i]<<std::endl;

output.close();

};//Get _N

int Vettore :: GetN () const {

return _N;

};

//Operatore di Assegnazione

Vettore & Vettore::operator =(const Vettore& vetty){

if (_v) delete [] _v;

_N=vetty.GetN();
_v=new double [_N];

for(int n; n<_N; n++)
_v[n]=vetty._v[n];

return *this;
};//Algoritmo Quicksort

void Vettore :: Sort (){
Quicksort(0,GetN()-1);
};

void Vettore :: Quicksort (unsigned int primo, unsigned int ultimo) {

if(ultimo-primo<=1){
if (GetComponent(primo)>GetComponent(ultimo)) Scambia(primo, ultimo);
return;
}

double pivot= GetComponent(int((primo+ultimo)/2));
unsigned int basso= primo, alto=ultimo;
while(basso < alto) {

while (GetComponent(basso)<pivot) basso++;
while (GetComponent(alto)>pivot) alto--;
if(basso<alto) { Scambia(basso,alto); basso++;};
};
Quicksort(primo, basso-1);
Quicksort(basso, ultimo);};

void Vettore :: Scambia(int a, int b){

double k;
k=_v[a];
_v[a]=_v[b];
_v[b]=k;
};

//Operatore somma fra vettori
VettoreLineare VettoreLineare::operator+ (const VettoreLineare& v){
VettoreLineare sum(v.GetN());

for(int i=0; i<_N; i++)
sum.SetComponent(i, _v[i]+v.GetComponent(i));

return sum;
};

//Operatore Moltiplicazione scalare

VettoreLineare VettoreLineare::operator* (double lambda){

for(int i=0; i<_N; i++)
_v[i]=_v[i]*lambda;

return *this;
};

//Operatore Assegnazione

VettoreLineare& VettoreLineare::operator= (const VettoreLineare& k){
if (_v) delete [] _v;
_N=k.GetN();
_v=new double [k.GetN()];

for (int i=0; i<_N; i++)
_v[i]=k.GetComponent(i);

return *this;
};

int VettoreLineare :: GetN() const {
return _N;
};

void VettoreLineare :: SetComponent(int i, double x) {

_v[i]=x;
};

main.c

#include "equazione_differenziale.h"#include "Vettore.h"
#include "iostream"
#include "TGraph.h"#include "TApplication.h"#include "TCanvas.h"#include "TAxis.h"
using namespace std;

int main () {//PRIMO PUNTO//Dichiarazione equazione
Protone* myProt=new Protone (1.67E-27, 1.60E-19, 1.E7, 0.2, 5.E8);
Runge_Kutta myKutta(myProt);

//Dichiarazione DatiIniziali
VettoreLineare DatiIniziali (2);
DatiIniziali.SetComponent(0, 0);
DatiIniziali.SetComponent(1,1E8);

//dichiarazione tempo
double t_min=0;
double h=1E-12;

//definizione variabili root
TApplication myApp("myApp",0,0);
TGraph* g = new TGraph();

//ciclo
for(int i=0;i<(1E-7-t_min)/h;i++){
double x,y;
x=t_min+i*h;
DatiIniziali= myKutta.Passo(x, DatiIniziali ,h);
y=DatiIniziali.GetComponent(0);
g->SetPoint(i,x,y);
};

//Run grafico
TCanvas *c=new TCanvas("C1","C1",1);
c->cd();
g->GetXaxis()->SetTitle("t[s]");
g->GetYaxis()->SetTitle("x[m]");
g->Draw("AL");
myApp.Run();return 0;

};

Странно то, что эта программа работает на университетских компьютерах, но не работает ни на одном из моих двух компьютеров. Я думаю, это означает, что я плохо установил ROOT на обоих компьютерах, но я искренне не знаю, как это доказать.

0

Решение

Вы заявляете

class EqDifferenzialeBase {

protected:
FunzioneBase* _f;

а также

class Runge_Kutta: public EqDifferenzialeBase {

protected:
FunzioneBase* _f;

а затем сделать

Runge_Kutta :: Runge_Kutta (FunzioneBase* f) : EqDifferenzialeBase(f) {};

а затем использовать его как

VettoreLineare Runge_Kutta::Passo(double t, VettoreLineare& v, double h){

VettoreLineare k1=_f->Eval(t,v);

Из этих строк должно быть довольно ясно, что идет не так и где segmentation violation берет свое начало. Так что избавьтесь от призрака, который затеняет поле, которое вы действительно хотите использовать.

И, как и в другом ответе, исправьте скалярное произведение векторного класса на стандарт оператора суммы.

1

Другие решения

Других решений пока нет …

По вопросам рекламы [email protected]