Runge-Kutta-4.0rdnung



  • Hallo Leute,

    ich habe ein kleines c++ code geschrieben, der das Runge-Kutta Verfahren anwendet. In dem Code muss man nur die Anfangsbedingungen und die zu lösende DGL eingeben. Ich habe vor kurzem mit dem Coden begonnen und bin offen für Verbesserungsvorschläge 🙂

    Mein Problem ist nun, dass wenn meine DGL nichtlinear wird, (z.B. y´=y^2) der Code nicht funktioniert (Die Werte gehen ins Unendliche). Bei der DGL (y´=y) funktioniert es sauber. Hier ist der Code:

    #include <iostream>
    #include <fstream>
    #include <math.h>
    #include <iomanip>
    #include <stdio.h>
    #include <stdlib.h>
    #include <string.h>
    
    
    using namespace std;
    
    double dgl(double x, double y);
    void rk();
    
    double ABx = 0;            // Anfangsbedingung für x
    double ABy = 1;            // Anfangsbedingung für y
    int anzahlPunkte = 50;     // Anzahl der Punkte, die geplottet werden 
    double h = 0.2;            // Schrittweite
    
    
    
    int main()
    {
    
    	dgl(ABx,ABy);
    	rk();
    
    return 0;
    }
    
    double dgl(double x, double y)
    {
    	double rslt;         // rslt == y´
    	rslt = y;
    	//rslt = y-2*x/y;
    	//rslt = (y*y);
    	return rslt;
    }
    
    void rk()
    {
    
    double k1,k2,k3,k4,m,x,y;
    double x_[anzahlPunkte], y_[anzahlPunkte];
    
    x = ABx;          // Anfangsbedingungen     
    y = ABy;
    
    	for (int i=1; i<anzahlPunkte; i++)
    	{
    	k1 = dgl(x, y);
    	k2 = dgl(x+h/2.0,(y+(h/2.0)*k1));
    	k3 = dgl(x+h/2.0,(y+(h/2.0)*k2));
    	k4 = dgl(x+h, (y+h*k3));
    	m = h/6*(k1 + 2*k2 + 2*k3 + k4);
    	x = x + h;
    	y = y + m; 
    
    	x_[i] = x;
    	y_[i] = y;
    	}
    
    		
    		for (int i=1; i<anzahlPunkte; i++)
    		{
    		cout << scientific << x_[i] << "  " << y_[i] << endl;
    		}
    		
    	ofstream file("rk.dat");
    	
    	file << scientific << ABx << " " << ABy << endl;
    	for (int i=1; i<anzahlPunkte; i++)
    	{
    	file << scientific << x_[i] << " " << y_[i] << endl;	
    	}
    	file.close();
    }
    


  • Würdest Du bitte Deinen Code nochmal ordentlich eingerückt in Deinen Beitrag einfügen und in eine Zeile vor Deinen Code ``` schreiben und in eine Zeile nach deinem Code auch? Alternativ kannst Du den eingefügten Code markieren, im Dropdown über dem Eingabefenster "C++" auswählen (sollte voreingestellt sein) und dann auf das </> daneben klicken.

    Den Menüpunkt "Bearbeiten" findest Du in dem Drei-Punkte-Menü rechts unter Deinem Beitrag.



  • @Swordfish Danke, hab den code richtig eingefügt. 🙂



  • @neoche dein Code ist eigentlich vollkommen korrekt. Dass inf herauskommt, liegt daran, dass ein double schlicht und ergreifend keine so großen Werte darstellen kann. long double kann (systemabhängig) ein wenig größer sein. Aber auch hier ist die Grenze bei y^2 schnell erreicht. Für größere Werte müsstest du dann spezielle Mathe-Libs verwenden (aber: brauchst du das wirklich?).

    Hier aber noch einmal dein Code mit ein paar kleinen Verschönerungen:

    #include <iostream>
    #include <fstream>
    #include <iomanip>
    
    static constexpr std::size_t amount_of_points = 50;
    static constexpr long double step_size = 0.2L;
    
    long double function(long double x, long double y)
    {
    //	return y - 2 * x / y;
    //	return y;
    	return y * y;
    }
    
    void output_point(std::ostream& os, std::pair<long double, long double> point)
    {
    	os << std::scientific << point.first << ' ' << point.second << '\n';
    }
    
    void rk(long double ABx, long double ABy)
    {
    	std::ofstream file("rk.dat");
    	output_point(std::cout, {ABx, ABy});
    	output_point(file, {ABx, ABy});
    
    	long double k1, k2, k3, k4, m;
    	for(std::size_t i = 1; i < amount_of_points; ++i)
    	{
    		k1 = function(ABx, ABy);
    		k2 = function(ABx + step_size/2.L, ABy + (step_size/2.L) * k1);
    		k3 = function(ABx + step_size/2.L, ABy + (step_size/2.L) * k2);
    		k4 = function(ABx + step_size, ABy + step_size * k3);
    		m = step_size/6.L * (k1 + 2 * k2 + 2 * k3 + k4);
    		ABx = ABx + step_size;
    		ABy = ABy + m;
    
    		output_point(std::cout, {ABx, ABy});
    		output_point(file, {ABx, ABy});
    	}
    }
    
    int main()
    {
    	rk(0, 1);
    	return 0;
    }
    
    


  • Das Verhalten dieser DGL ist eben so, da nützt es auch nichts, größere Datentypen zu nehmen.



  • @Unterfliege Ich danke Dir erstmals!! Dann liegt es doch nicht am code, sondern an der von mir hergeleitete Gleichung! Außerdem hat die Gleichung eine Singularität bei x=1, wenn ich mich nicht irre. Danke für die Verbesserungsvorschläge 🙂



  • @Bashar Da hast du wohl Recht!!


Anmelden zum Antworten