Trapez-Integration, Romberg-Verfahren
-
Hallo!
Eine Aufgabe meines C++ Kurses, Trapez-Integration, Romberg-Verfahren..
..eigentlich habe ich die Aufgabe gelöst, ich habe nur ein kleines Problem, die genauigkeit!
folgender Code :
#include <iostream> #include <fstream> #include <iomanip> #include <cmath> #include <string> #include <vector> #include <cstdlib> #include <new> using namespace std; //_________________________________Funktionen für Integration___________________________________// //Integrier Klasse class integ { //__________________________________________________Private private: double a; double b; //Übernahme der Funktion FUNKTIONIERT! void (*P)(double&); double f(double x,void (*pt2Func)(double&)) {pt2Func(x);return(x);} //Funktionswert an der Stelle j FUNKTIONIERT! double F_j (double j,double N) {return(f(a+j*(b-a)/N,P));} //Trapezsumme FUNKTIONIERT! double Ftr_N (double N) { double Sum=F_j(0,N)+F_j(N,N); for (int i=1;i<N;i++) Sum+=(2*F_j(i,N)); Sum *= (b-a)/(N*2); //cout << "\n N:" << N << " Sum" << Sum; return(Sum); } //Interpolation double ipol(double Pio[] ,int Len=0) { cout << "\n\n"; long FOR=4; //Ausgabe der ersten Zeile for (int j=Len-1;j>=0;j--) cout << setw(3) << "[" << j << ".0]" << setprecision(4) << setw(5) << Pio[j]; cout << "\n"; //Interpolierter Wert for (int h=1;h<Len;h++) { for (int j=Len-1;j>=h;j--) { Pio[j]=Pio[j]+(Pio[j]-Pio[j-1])/(FOR-1); cout << setw(3) << "[" << j << "." << h << "]" << setprecision(3) << setw(5) << Pio[j]; } FOR*=4; cout << "\n"; } return(Pio[Len-1]); //return(0); } //Eigentliche Integration double intgr() { const double e=1E-4; int i=2; double T=0; double G=Ftr_N(1); //ermittle N oder auch i max for (;abs((G-T)/G)>e and i<10/*65536*/;i*=2) { T=G; G=Ftr_N(i); } //erstelle die Pio in einem Array i++; //i intervalle, aber i+1 Stützstellen double Pio[i]; for (int j=1;j<i;j++) { Pio[j]=Ftr_N(2*j); } Pio[0]=0;//dieser Wert ändert NIX! ..ob 100 oder 200.. völlig egal //Romberg interpolation double xpol = ipol(Pio,i); return (xpol); } //__________________________________________________Public public: //Konstruktor integ(void (*pt2Func)(double&),double A=0,double B=0) { P = *pt2Func; a=A; b=B; } //Ausgabe void Solve() { cout << "\n Integration von f auf dem Intervall [ " << a << " ; " << b << " ] ergibt S = " << setprecision(10) << intgr(); } }; //_________________________________Main___________________________________// //Die Funktion über die integriert wird void f(double &x) {x=sqrt(x)/(1+x*x);} //void f(double &x) {x=x*x;} //Testfunktion int main() { int a=0,b=10; cout << "\n 1) Trapez-Integration, Romberg-Verfahren : \n"; integ INT(f,a,b); INT.Solve(); cout << "\n\n"; return(0); }Die verwendete Funktion f(x) = sqrt(x)/(1+x²) solte von 0 bis 10 integriert ungefähr 1,590243869 ergeben! ich kriege mit meinem Code aber nur 1.551838253!
..ich habe absolut keinen Plan wo der Fehler liegt, vielleicht ist das Verfahren auch nicht so präzise, aber schließlich hat man uns ja so viele Stellen nach dem Komma gegeben..
Ich hab schon ziemlich alles nachgerechnet, teilweise per Hand, und keinen Fehler gefunden
HILFE!! x)
danke für eure Zeit
Victor
-
ich hab' jetzt keine lust, das ganze auseinanderzufrickeln. aber zwei tipps:
- schau dir, wie der wert mit steigendem N konvergiert (oder auch nicht)
also N = 2,20,200,2000... und dann differenz zum exakten wert plotten. bei deinem integrationsverfahren sollte ja eine konvergenzgeschwindigkeit angegeben sein
- löse ein paar trivial integrierbare probleme (f(x)=1, x^n, exp(x)) u.s.w. dann siehst du vielleicht den fehler.//Die Funktion über die integriert wird void f(double &x) {x=sqrt(x)/(1+x*x);}warum verwendest du nicht
//Die Funktion über die integriert wird double f(double x) {return sqrt(x)/(1+x*x);}dann kannst du dir auch integ::f(...) sparen?
-
EDIT 2 : gelößt! (kann man den thread damit markieren?) die FOR variable muss als datantyp long double haben! dann gehts
EDIT : liegt wohl an der 4^k in der Extrapolation, die wird sehr schnell viel zu groß und wird somit 0.. muss die extrapolation anders berechnen.
Danke für deine Antwort,
das mit der Funktion hab ich geändert, hatte es anders gemacht am anfang da F_j jetzt einen Parameter mehr hat aber egal, so spare ich immerhin eine Funktion _
#include <iostream> #include <fstream> #include <iomanip> #include <cmath> #include <string> #include <vector> #include <cstdlib> #include <new> using namespace std; //_________________________________Funktionen für Integration___________________________________// //Integrier Klasse class integ { //__________________________________________________Private private: double a; double b; //Übernahme der Funktion FUNKTIONIERT! double (*P)(double); //Funktionswert an der Stelle j FUNKTIONIERT! double F_j (double j,double N,double (*pt2Func)(double)) {return(pt2Func(a+j*(b-a)/N));} //Trapezsumme FUNKTIONIERT! double Ftr_N (double N) { double Sum=F_j(0,N,P)+F_j(N,N,P); for (int i=1;i<N;i++) Sum+=(2*F_j(i,N,P)); Sum *= (b-a)/(N*2); cout << "\n N:" << N << " Sum" << Sum; return(Sum); } //Interpolation double ipol(double Pio[] ,int Len=0) { cout << "\n\n"; long FOR=4; //Ausgabe der ersten Zeile //for (int j=Len-1;j>=0;j--) cout << setw(3) << "[" << j << ".0]" << setprecision(4) << setw(5) << Pio[j]; cout << "\n"; //Interpolierter Wert for (int h=1;h<Len;h++) { for (int j=Len-1;j>=h;j--) { Pio[j]=Pio[j]+(Pio[j]-Pio[j-1])/(FOR-1); //cout << setw(3) << "[" << j << "." << h << "]" << setprecision(3) << setw(5) << Pio[j]; } FOR*=4; //cout << "\n"; } return(Pio[Len-1]); //return(0); } //Eigentliche Integration double intgr() { const double e=1E-4; int i=2; double T=0; double G=Ftr_N(1); //ermittle N oder auch i max for (;abs((G-T)/G)>e and i<1000/*65536*/;i*=2) { T=G; G=Ftr_N(i); } //erstelle die Pio in einem Array i++; //i intervalle, aber i+1 Stützstellen double Pio[i]; for (int j=1;j<i;j++) { Pio[j]=Ftr_N(2*j); } Pio[0]=0;//dieser Wert ändert NIX! ..ob 100 oder 200.. völlig egal //Romberg interpolation double xpol = ipol(Pio,i); return (xpol); } //__________________________________________________Public public: //Konstruktor integ(double (*pt2Func)(double),double A=0,double B=0) { P = *pt2Func; a=A; b=B; } //Ausgabe void Solve() { cout << "\n Integration von f auf dem Intervall [ " << a << " ; " << b << " ] ergibt S = " << setprecision(10) << intgr(); } }; //_________________________________Main___________________________________// //Die Funktion über die integriert wird double f(double x) {return sqrt(x)/(1+x*x);} //void f(double &x) {x=x*x;} //Testfunktion int main() { int a=0,b=10; cout << "\n 1) Trapez-Integration, Romberg-Verfahren : \n"; integ INT(f,a,b); INT.Solve(); cout << "\n\n"; return(0); }Also die Trapezsummen konvergieren für N -> ∞ gegen den richtigen Wert, das funktioniert anscheinend sehr genau, aber den Wert nehme ich ja nur als erste Zeile im Tableau der Extrapolation.
Der per Extrapolation berechnete Wert konvergiert jedoch gegen den Falschen Wert
Mit einfachen Funktion (1 x x² ..) merkt man davon nichts, erst bei Funktion mit mehreren Termen wie x²+3 +x
So verhält sich die letzte Zeile meines tableaus für ein relativ kleinen N, der letzte Wert müsste der gesuchten Integration entsprechen.
[512.6] 1.59
[512.7] 1.59
[512.8] 1.59
[512.9] 1.59
[512.10] 1.59
[512.11] 1.59
[512.12] 1.59
[512.13] 1.59
[512.14] 1.59
[512.15] 1.59
[512.16] 1.59
[512.17] 1.59
[512.18] 1.59
[512.19] 1.59
[512.20] 1.59
[512.21] 1.59
[512.22] 1.59
[512.23] 1.59
[512.24] 1.59
[512.25] 1.59
[512.26] 1.59
[512.27] 1.59
[512.28] 1.59
[512.29] 1.59
[512.30] 1.59
[512.31] 1.59
[512.32] 1.59
[512.33] 1.59
[512.34] 1.59
[512.35] 1.59
[512.36] 1.59
[512.37] 1.59
[512.38] 1.59
[512.39] 1.59
[512.40] 1.59
[512.41] 1.59
[512.42] 1.59
[512.43] 1.59
[512.44] 1.59
[512.45] 1.59
[512.46] 1.59
[512.47] 1.59
[512.48] 1.59
[512.49] 1.59
[512.50] 1.59
[512.51] 1.59
[512.52] 1.59
[512.53] 1.59
[512.54] 1.59
[512.55] 1.59
[512.56] 1.59
[512.57] 1.59
[512.58] 1.59
[512.59] 1.59
[512.60] 1.59
[512.61] 1.59
[512.62] 1.59
[512.63] 1.59
[512.64] 1.59
[512.65] 1.59
[512.66] 1.59
[512.67] 1.59
[512.68] 1.59
[512.69] 1.59
[512.70] 1.59
[512.71] 1.59
[512.72] 1.59
[512.73] 1.59
[512.74] 1.59
[512.75] 1.59
[512.76] 1.59
[512.77] 1.59
[512.78] 1.59
[512.79] 1.59
[512.80] 1.59
[512.81] 1.59
[512.82] 1.59
[512.83] 1.59
[512.84] 1.59
[512.85] 1.59
[512.86] 1.59
[512.87] 1.59
[512.88] 1.59
[512.89] 1.59
[512.90] 1.59
[512.91] 1.59
[512.92] 1.59
[512.93] 1.59
[512.94] 1.59
[512.95] 1.59
[512.96] 1.59
[512.97] 1.59
[512.98] 1.59
[512.99] 1.59
[512.100] 1.59
[512.101] 1.59
[512.102] 1.59
[512.103] 1.59
[512.104] 1.59
[512.105] 1.59
[512.106] 1.59
[512.107] 1.59
[512.108] 1.59
[512.109] 1.59
[512.110] 1.59
[512.111] 1.59
[512.112] 1.59
[512.113] 1.59
[512.114] 1.59
[512.115] 1.59
[512.116] 1.59
[512.117] 1.59
[512.118] 1.59
[512.119] 1.59
[512.120] 1.59
[512.121] 1.59
[512.122] 1.59
[512.123] 1.59
[512.124] 1.59
[512.125] 1.59
[512.126] 1.59
[512.127] 1.59
[512.128] 1.59
[512.129] 1.59
[512.130] 1.59
[512.131]1.5899
[512.132]1.5899
[512.133]1.5899
[512.134]1.5899
[512.135]1.5899
[512.136]1.5899
[512.137]1.5899
[512.138]1.5899
[512.139]1.5899
[512.140]1.5899
[512.141]1.5899
[512.142]1.5899
[512.143]1.5899
[512.144]1.5899
[512.145]1.5899
[512.146]1.5899
[512.147]1.5899
[512.148]1.5899
[512.149]1.5899
[512.150]1.5899
[512.151]1.5899
[512.152]1.5899
[512.153]1.5899
[512.154]1.5899
[512.155]1.5899
[512.156]1.5899
[512.157]1.5899
[512.158]1.5899
[512.159]1.5899
[512.160]1.5899
[512.161]1.5899
[512.162]1.5899
[512.163]1.5899
[512.164]1.5899
[512.165]1.5899
[512.166]1.5899
[512.167]1.5899
[512.168]1.5899
[512.169]1.5899
[512.170]1.5899
[512.171]1.5899
[512.172]1.5899
[512.173]1.5899
[512.174]1.5899
[512.175]1.5899
[512.176]1.5899
[512.177]1.5899
[512.178]1.5899
[512.179]1.5899
[512.180]1.5899
[512.181]1.5899
[512.182]1.5899
[512.183]1.5899
[512.184]1.5899
[512.185]1.5899
[512.186]1.5899
[512.187]1.5899
[512.188]1.5899
[512.189]1.5899
[512.190]1.5899
[512.191]1.5899
[512.192]1.5899
[512.193]1.5899
[512.194]1.5899
[512.195]1.5899
[512.196]1.5899
[512.197]1.5899
[512.198]1.5899
[512.199]1.5899
[512.200]1.5899
[512.201]1.5898
[512.202]1.5898
[512.203]1.5898
[512.204]1.5898
[512.205]1.5898
[512.206]1.5898
[512.207]1.5898
[512.208]1.5898
[512.209]1.5898
[512.210]1.5898
[512.211]1.5898
[512.212]1.5898
[512.213]1.5898
[512.214]1.5898
[512.215]1.5898
[512.216]1.5898
[512.217]1.5898
[512.218]1.5898
[512.219]1.5898
[512.220]1.5898
[512.221]1.5898
[512.222]1.5898
[512.223]1.5898
[512.224]1.5898
[512.225]1.5898
[512.226]1.5898
[512.227]1.5898
[512.228]1.5898
[512.229]1.5898
[512.230]1.5898
[512.231]1.5898
[512.232]1.5898
[512.233]1.5898
[512.234]1.5898
[512.235]1.5898
[512.236]1.5898
[512.237]1.5898
[512.238]1.5898
[512.239]1.5898
[512.240]1.5898
[512.241]1.5898
[512.242]1.5898
[512.243]1.5898
[512.244]1.5898
[512.245]1.5898
[512.246]1.5898
[512.247]1.5897
[512.248]1.5897
[512.249]1.5897
[512.250]1.5897
[512.251]1.5897
[512.252]1.5897
[512.253]1.5897
[512.254]1.5897
[512.255]1.5897
[512.256]1.5897
[512.257]1.5897
[512.258]1.5897
[512.259]1.5897
[512.260]1.5897
[512.261]1.5897
[512.262]1.5897
[512.263]1.5897
[512.264]1.5897
[512.265]1.5897
[512.266]1.5897
[512.267]1.5897
[512.268]1.5897
[512.269]1.5897
[512.270]1.5897
[512.271]1.5897
[512.272]1.5897
[512.273]1.5897
[512.274]1.5897
[512.275]1.5897
[512.276]1.5897
[512.277]1.5897
[512.278]1.5897
[512.279]1.5897
[512.280]1.5896
[512.281]1.5896
[512.282]1.5896
[512.283]1.5896
[512.284]1.5896
[512.285]1.5896
[512.286]1.5896
[512.287]1.5896
[512.288]1.5896
[512.289]1.5896
[512.290]1.5896
[512.291]1.5896
[512.292]1.5896
[512.293]1.5896
[512.294]1.5896
[512.295]1.5896
[512.296]1.5896
[512.297]1.5896
[512.298]1.5896
[512.299]1.5896
[512.300]1.5896
[512.301]1.5896
[512.302]1.5896
[512.303]1.5896
[512.304]1.5895
[512.305]1.5895
[512.306]1.5895
[512.307]1.5895
[512.308]1.5895
[512.309]1.5895
[512.310]1.5895
[512.311]1.5895
[512.312]1.5895
[512.313]1.5895
[512.314]1.5895
[512.315]1.5895
[512.316]1.5895
[512.317]1.5895
[512.318]1.5895
[512.319]1.5895
[512.320]1.5895
[512.321]1.5895
[512.322]1.5895
[512.323]1.5894
[512.324]1.5894
[512.325]1.5894
[512.326]1.5894
[512.327]1.5894
[512.328]1.5894
[512.329]1.5894
[512.330]1.5894
[512.331]1.5894
[512.332]1.5894
[512.333]1.5894
[512.334]1.5894
[512.335]1.5894
[512.336]1.5894
[512.337]1.5894
[512.338]1.5894
[512.339]1.5893
[512.340]1.5893
[512.341]1.5893
[512.342]1.5893
[512.343]1.5893
[512.344]1.5893
[512.345]1.5893
[512.346]1.5893
[512.347]1.5893
[512.348]1.5893
[512.349]1.5893
[512.350]1.5893
[512.351]1.5893
[512.352]1.5892
[512.353]1.5892
[512.354]1.5892
[512.355]1.5892
[512.356]1.5892
[512.357]1.5892
[512.358]1.5892
[512.359]1.5892
[512.360]1.5892
[512.361]1.5892
[512.362]1.5892
[512.363]1.5891
[512.364]1.5891
[512.365]1.5891
[512.366]1.5891
[512.367]1.5891
[512.368]1.5891
[512.369]1.5891
[512.370]1.5891
[512.371]1.5891
[512.372]1.589
[512.373]1.589
[512.374]1.589
[512.375]1.589
[512.376]1.589
[512.377]1.589
[512.378]1.589
[512.379]1.589
[512.380]1.5889
[512.381]1.5889
[512.382]1.5889
[512.383]1.5889
[512.384]1.5889
[512.385]1.5889
[512.386]1.5889
[512.387]1.5888
[512.388]1.5888
[512.389]1.5888
[512.390]1.5888
[512.391]1.5888
[512.392]1.5888
[512.393]1.5888
[512.394]1.5887
[512.395]1.5887
[512.396]1.5887
[512.397]1.5887
[512.398]1.5887
[512.399]1.5886
[512.400]1.5886
[512.401]1.5886
[512.402]1.5886
[512.403]1.5886
[512.404]1.5885
[512.405]1.5885
[512.406]1.5885
[512.407]1.5885
[512.408]1.5885
[512.409]1.5884
[512.410]1.5884
[512.411]1.5884
[512.412]1.5884
[512.413]1.5883
[512.414]1.5883
[512.415]1.5883
[512.416]1.5883
[512.417]1.5882
[512.418]1.5882
[512.419]1.5882
[512.420]1.5882
[512.421]1.5881
[512.422]1.5881
[512.423]1.5881
[512.424]1.588
[512.425]1.588
[512.426]1.588
[512.427]1.5879
[512.428]1.5879
[512.429]1.5879
[512.430]1.5878
[512.431]1.5878
[512.432]1.5877
[512.433]1.5877
[512.434]1.5877
[512.435]1.5876
[512.436]1.5876
[512.437]1.5875
[512.438]1.5875
[512.439]1.5874
[512.440]1.5874
[512.441]1.5873
[512.442]1.5873
[512.443]1.5872
[512.444]1.5872
[512.445]1.5871
[512.446]1.5871
[512.447]1.587
[512.448]1.587
[512.449]1.5869
[512.450]1.5868
[512.451]1.5868
[512.452]1.5867
[512.453]1.5866
[512.454]1.5865
[512.455]1.5865
[512.456]1.5864
[512.457]1.5863
[512.458]1.5862
[512.459]1.5861
[512.460]1.586
[512.461]1.5859
[512.462]1.5858
[512.463]1.5857
[512.464]1.5856
[512.465]1.5855
[512.466]1.5854
[512.467]1.5853
[512.468]1.5852
[512.469]1.585
[512.470]1.5849
[512.471]1.5848
[512.472]1.5846
[512.473]1.5845
[512.474]1.5843
[512.475]1.5841
[512.476]1.5839
[512.477]1.5837
[512.478]1.5836
[512.479]1.5833
[512.480]1.5831
[512.481]1.5829
[512.482]1.5826
[512.483]1.5824
[512.484]1.5821
[512.485]1.5818
[512.486]1.5815
[512.487]1.5812
[512.488]1.5808[512.489]1.5805
[512.490]1.5801
[512.491]1.5797
[512.492]1.5792
[512.493]1.5787
[512.494]1.5782
[512.495]1.5776
[512.496]1.577
[512.497]1.5764
[512.498]1.5757
[512.499]1.5749
[512.500]1.574
[512.501]1.5731
[512.502]1.5721
[512.503]1.571
[512.504]1.5697
[512.505]1.5683
[512.506]1.5668
[512.507]1.565
[512.508]1.5631
[512.509]1.5608
[512.510]1.5583
[512.511]1.5553
[512.512]1.5518Integration von f auf dem Intervall [ 0 ; 10 ] ergibt S = 1.551838253
gibt es eine möglichkeit so eine lange liste so zu posten das man das auf und zuklappen kann? ..so wie etwa ein spoiler ..nur um platz zu sparen