Doppelpendel-DGL ohne Reibung und ohne Näherungen.
Befestigungspunkt - Stange l1 - Masse m1 - Stange l2 - Masse m2.
Beschreibung durch die beiden Winkel der Stangen (phi, psi), Lagrange-Ansatz, etc.
In den Bildren aufgetragen ist die Lage der ("unteren") Masse m2 in der Ebene; die Schwerkraft wirkt nach unten.
Anfangsbedingung: Alle Simulationen (ausser anders erwähnt) wurden mit den Anfangsbedingungen phi=0.9*Pi, psi=Pi oder mit phi=0.9*Pi, psi=0.5*Pi durchgeführt (welche es sind, kann man unschwer an der Lage des Anfangs der Trajektorie erkennen).

Mittels adaptivem explizitem Einschrittverfahren der Ordnung 2 (Heun, Collatz) und 4 (Runge-Kutta).
Vergleich der verschiedenen Verfahren bei adaptiver Schrittweitenkontrolle und einem maximalen Diskretisierungsfehler eps=1e-10. (Die Verfahren niedrigerer Ordnung brauchen dazu natürlich kleinere Schrittweiten und dementsprechend länger.)
Man beachte rechts unten: Es handelt sich um ein und das selbe Verfahren aber zwei verschiedenen Implementierungen.

Vergleich des Runge-Kutta-Verfahrens ("klassisches", also ältestes RK-Verfahren) bei verschiedenen Diskretisierungsfehler-Schranken eps=1e-5 bis eps=1e-13. Wie immer adaptives Verfahren mit gleichen Anfangsbedingungen.

Die eine Trajektorie benutzt die selben Anfangsbedingungen wie oben, bei der anderen ist der Start-Winkel von m1 von 0.9*Pi auf 0.9*Pi+1e-12 gestört (also eine relative Änderung um etwa 3e-13).
Simulation mit verschiedenen Fehlerschranken eps.

Rechts: "Klassisches" Runge-Kutta-Verfahren mit Diskretisierungsfehler eps=1e-9; verwendete Schrittweite für Doppelschritt unten über die Zeit aufgetragen. (Das verwendete adaptive Verfahren berechnet zwei Schritte mit halber Schrittweite (="Doppelschritt") und einen mit voller Schrittweite und schätzt den Diskretisierungsfehler aus der Differenz der Ergebnisse. Ist diese klein genug, so wird das Ergebnis des Doppelschrittes akzeptiert.)

This code is Copyright (c) 02/2004 by Wolfgang Wieser.
It may be copied under the terms of the GNU GPL.

The complete ODE code is part of AniVision, version 0.3.4a or above.
The graphical test program for the ODE solver routines is available here: odetest.cc.
[Requires QTXlib (bundeled with burn simulation), hlib, Qt and probably Linux and GCC -- and AniVision, of course]

The ODE function:

class MyODEFunc : public NUM::Function_ODE
{
	protected:
		// [overriding virtual]
		void function(double /*x*/,const double *y,double *result)
		{
			register double phi=y[0],phip=y[1],psi=y[2],psip=y[3];
			
			double tmp = m2*l1*l2*sin(phi-psi);
			double a = -tmp*psip*psip - g*(m1+m2)*l1*sin(phi);
			double d =  tmp*phip*phip - g*    m2 *l2*sin(psi);
			double b = l1*l1*(m1+m2);
			double c = m2*l1*l2*cos(phi-psi);
			double e = l2*l2*m2;
			double det = (b*e-c*c);
			
			result[0]=phip;
			result[1]=(a*e-c*d)/det;
			result[2]=psip;
			result[3]=(b*d-c*a)/det;
			
		 #if 0
			// Control: calc energy:
			double T = 0.5*(m1+m2)*SQR(l1*phip) +
				m2*( 0.5*SQR(l2*psip) + l1*l2*phip*psip*cos(phi-psi) );
			double U = -g*( (m1+m2)*l1*cos(phi) + m2*l2*cos(psi) );
			double E = T+U;
			
			static int cnt=0;
			if(!(++cnt%100))
			{  fprintf(stderr,"\n%g",E);  }
		 #endif
		}
	public:
		MyODEFunc() : Function_ODE(4) {}
		~MyODEFunc() {}
};

The solver algorithm (methods below):

void ODESolver::SetAdaptivePars(double epsilon)
{
	// Get order of consistency:
	int p=GetSolverParams()->p;
	
	adapt_limit1 = (1.0-pow(2.0,-p)) * epsilon;
	adapt_limit0 = adapt_limit1 * pow(2.0,p+1);
	// Add some safety margin to prevent up-down-up-down-up - switches.
	// This should normally not be necessary.
	adapt_limit1 *= 0.9;
}


void ODESolver::SingleStep(Function_ODE &odefunc,
	double &xn,double h,double *yn)
{
	ComputeStep(odefunc,xn,h,yn);
	++nsteps;
	xn+=h;
	if(used_h_min>h)  used_h_min=h;
	if(used_h_max<h)  used_h_max=h;
}


void ODESolver::AdaptiveStep(Function_ODE &odefunc,
	double &xn,double &h,double *yn)
{
	int dim=odefunc.Dim();
	double ynF[dim],ynH[dim],ynHH[dim];
	
	double ah=fabs(h);
	for(bool first=1,last=(ah<=h_min);;first=0)
	{
		if(first)
		{
			for(int d=0; d<dim; d++)
			{  ynF[d]=yn[d];  ynH[d]=yn[d];  }
			
			// Do the full step:
			ComputeStep(odefunc,xn,h,ynF);
		}
		else
		{
			// Use saved half step from last time as full step:
			for(int d=0; d<dim; d++)
			{  ynF[d]=ynHH[d];  ynH[d]=yn[d];  }
		}
		
		if(used_h_min>ah)  used_h_min=ah;
		if(used_h_max<ah)  used_h_max=ah;
		
		// Do the two half steps:
		double h2=0.5*h;
		ComputeStep(odefunc,xn,h2,ynH);
		
		// Save that in case we need to reduce the step size:
		if(!last)  for(int d=0; d<dim; d++)  ynHH[d]=ynH[d];
		
		ComputeStep(odefunc,xn+h2,h2,ynH);
		
		// Compute the difference:
		// (maximum norm of full_step minus two_half_steps)
		double delta=0.0;
		for(int d=0; d<dim; d++)
		{
			double diff=fabs(ynF[d]-ynH[d]);
			if(!d || delta<diff)  delta=diff;
		}
		
		if(delta<=adapt_limit0)
		{
			// Accept result.
			xn+=h;
			
			if(delta<=adapt_limit1)
			{
				// Even try larger step size next time:
				h+=h;
				if(h>0.0)
				{  if(h>h_max)  h=h_max;  }
				else
				{  if(-h>h_max)  h=-h_max;  }
			}
			
			break;
		}
		
		// Must reduce step size.
		if(last)
		{  // Already smallest step size.
			xn+=h;
			break;
		}
		
		h=h2;
		if(h>0.0)
		{
			if(h<h_min)
			{  h=h_min;  last=1;  }
		}
		else
		{
			if(-h<h_min)
			{  h=-h_min;  last=1;  }
		}
	}
	
	nsteps+=2;
	for(int d=0; d<dim; d++)
		yn[d]=ynH[d];
}

The different methods:

void ODESolver_Heun::ComputeStep(Function_ODE &odefunc,
	double xn,double h,double *yn)
{
	int dim=odefunc.Dim();
	double k[2][dim];
	double ytmp[dim];
	
	odefunc(xn,yn,k[0]);
	
	for(int d=0; d<dim; d++)
		ytmp[d]=yn[d]+h*k[0][d];
	
	odefunc(xn+h,ytmp,k[1]);
	
	double h2=0.5*h;
	for(int d=0; d<dim; d++)
		yn[d]+=h2*(k[0][d]+k[1][d]);
}

void ODESolver_Collatz::ComputeStep(Function_ODE &odefunc,
	double xn,double h,double *yn)
{
	int dim=odefunc.Dim();
	double k[dim];
	double ytmp[dim];
	
	odefunc(xn,yn,k);
	
	double h2=0.5*h;
	for(int d=0; d<dim; d++)
		ytmp[d]=yn[d]+h2*k[d];
	
	odefunc(xn+h2,ytmp,k);
	
	for(int d=0; d<dim; d++)
		yn[d]+=h*k[d];
}

void ODESolver_RungeKutta::ComputeStep(Function_ODE &odefunc,
	double xn,double h,double *yn)
{
	int dim=odefunc.Dim();
	double k[dim];
	double ydelta[dim];
	double ytmp[dim];
	
	odefunc(xn,yn,k);
	
	double h2=0.5*h;
	for(int d=0; d<dim; d++)
	{
		ydelta[d]=k[d];
		ytmp[d]=yn[d]+h2*k[d];
	}
	
	odefunc(xn+h2,ytmp,k);
	
	for(int d=0; d<dim; d++)
	{
		ydelta[d]+=2.0*k[d];
		ytmp[d]=yn[d]+h2*k[d];
	}
	
	odefunc(xn+h2,ytmp,k);
	
	for(int d=0; d<dim; d++)
	{
		ydelta[d]+=2.0*k[d];
		ytmp[d]=yn[d]+h*k[d];
	}
	
	odefunc(xn+h,ytmp,k);
	
	h2=h/6.0;
	for(int d=0; d<dim; d++)
		yn[d]+=h2*(ydelta[d]+k[d]);
}

void ODESolver_GeneralRungeKutta::ComputeStep(Function_ODE &odefunc,
	double xn,double h,double *yn)
{
	int dim=odefunc.Dim();
	double k[sp.m][dim];
	double ytmp[dim];
	
	// First evaluation (optimized):
	odefunc(xn,yn,k[0]);
	
	// All the other evaluations:
	for(int i=1; i<sp.m; i++)
	{
		for(int d=0; d<dim; d++)
		{
			double tmp=0.0;
			for(int s=0; s<i; s++)
				tmp+=beta(i,s)*k[s][d];
			ytmp[d]=yn[d]+h*tmp;
		}
		
		odefunc(xn+alpha[i]*h,ytmp,k[i]);
	}
	
	// Update values:
	for(int d=0; d<dim; d++)
	{
		double delta=0.0;
		for(int s=0; s<sp.m; s++)
			delta+=gamma[s]*k[s][d];
		yn[d]+=h*delta;
	}
}