c++ - c++中的微分方程，Runge Kutta方法的问题

Question

我用 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;
};

主.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，但我真的不知道如何证明它。

Answer 1

你声明

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 的来源。因此，摆脱影响您真正想要使用的领域的幽灵。

和另一个答案一样，将 vector 类的标量积更正为求和运算符的标准。