C ++中的微分方程，Runge Kutta方法的问题 [英] differential equations in c++, problems with Runge Kutta's Method
本文介绍了C ++中的微分方程，Runge Kutta方法的问题的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我用c ++编写了一个程序，可以解决微分方程.问题是，它似乎不能与ROOT一起很好地工作.它可以编译，但是当我执行时，这就是我得到的:
I wrote a program in c++ that should solve differential equations. The problem is, it seems like it does not work well with ROOT. It compilates fine, but when I execute, this is what I get:
*** 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
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
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
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
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，但是我真诚地不知道该如何解决.
The strange thing is that this program works on university computers, but it doesn't on neither of my two computers. I'm thinking this means I installed ROOT badly on both computers, but I sincerely wouldn't know how to proove it.
推荐答案
您声明
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违规的起源.因此，请消除阴影，它会遮盖您真正想要使用的字段.
From these lines it should be pretty clear what goes wrong and where the segmentation violation originates. So get rid of the ghost that shadows the field you really want to use.
就像在其他答案中一样，将向量类的标量积校正为和运算符的标准.
And as in the other answer, correct the scalar product of the vector class to the standard of the sum operator.
这篇关于C ++中的微分方程，Runge Kutta方法的问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
