Adaptive step size algorithm in c++ for runge kutta

  • C/C++
  • Thread starter FunkyDwarf
  • Start date
  • #1
FunkyDwarf
489
0
Hey guys,

Attempting to write an adaptive step size function into a 4th order runge kutta integrator for basic orbits.

The problems (so far) are as follows:
1.) The step size does not seem to change as one would expect if the recursive definition of the StepSize function was working properly

2.) When compared with the mathematica solution, it would seem the orbit doesn't extend as far as it should in time, ie the arc length is shorter. It would seem this is due to the coordinates remaining constant after a certain point.

Code can be found here:
http://members.iinet.net.au/~housewrk/programtest.cpp [Broken]

Apologies in advance for stupid errors and crappy coding :)

Cheers
-G
 
Last edited by a moderator:

Answers and Replies

  • #2
Sagekilla
19
0
It looks like you're integrating slightly wrong. I wouldn't separate the Runge Kutta integration into x and y steps. You should solve for k1x = f(t, x) and k1y = f(t, y). Then, calculate k2x and k2y, and so forth.

The idea is you want to make a complete approximation going forward. You're stepping the components forward separately, but they're dependent on each other so you're introducing errors. For small enough stepsizes it might not matter, but it's still a bad idea.

Otherwise, for what you're doing, using a full step vs two half steps for error control is reasonable but I'd recommend taking a look into embedded Runge Kutta methods as they have a very easy way to implement error control.

I wrote up a ode integrator before and it's fairly easy to do this:
Code:
// Bogacki-Shampine method
void integrate(double &t, double &y, double &h)
{
	double k1, k2, k3, k4, y1, y2;

	// f(t, y) is our derivative
	k1 = f(t, y)
	k2 = f(t + h * 0.50, y + h * 0.50 * k1);
	k3 = f(t + h * 0.75, y + h * 0.75 * k2);

	// 3rd order accurate solution
	y1 = y + h * (2.0/9.0 * k1 + 1.0/3.0 * k2 + 4.0/9.0 * k3);

	k4 = f(t + h, y1);

	// 2nd order accurate solution
	y2 = y + h1 * (7.0/24.0 * k1 + 1.0/4.0 * k2 + 1.0/3.0 * k3 + 1.0/8.0 * k4);

	y = y2;
	t += h;

	// Define atol (absolute tolerance) and rtol (relative tolerance) somewhere
	// We use the difference between the two solutions as an approximation of the error, scaled by our desired absolute and relative tolerance
	// Both atol and rtol should be non-zero.

	double scale = atol + max(fabs(y1), fabs(y2)) * rtol;
	double error = max(fabs(y2 - y1), 2E-10) / scale;
	h = h * pow(error, -1.0/3.0);
}

Using something like an embedded method may be easier to do. You'll still want to control the integration and make sure a step doesn't introduce too big of an error (which the above code doesn't account for).

I'll have to look through your code some more and pick it apart to figure out what's wrong, though.
 
  • #3
rcgldr
Homework Helper
8,791
580
Even if you use an iterative method, such as predictor-corrector:

http://en.wikipedia.org/wiki/Predictor-corrector_method#Euler_Trapezoidal_Example

The algorithm may converge to a specific value for each step, but in that wiki example, it's a linear approximation for each step, and the overall error will increase with step size, regardless of convergence at each step. The method you show splits up each step into 4 steps, but it's still a linear weighted average over those 4 steps. You could just reduce step size by a factor of 4 and use a iterative method for more accuracy, but it would take more time. On a modern PC, I doubt CPU time is an issue for most projects.

If you want to vary the step size, it should be some inverse function of the magnitude of the function(s) you're trying to integrate, but this could lead to errors if the magnitude increased rapidly within a step period.

Wiki includes a good example of a "stiff equation" you can use to test your method with, since it's directly solvable:

http://en.wikipedia.org/wiki/Stiff_equation
 
Last edited:

Suggested for: Adaptive step size algorithm in c++ for runge kutta

Replies
15
Views
998
  • Last Post
Replies
0
Views
294
Replies
2
Views
442
  • Last Post
Replies
8
Views
1K
Replies
6
Views
507
Replies
0
Views
372
  • Last Post
Replies
13
Views
1K
Replies
10
Views
1K
  • Last Post
Replies
0
Views
311
  • Last Post
Replies
14
Views
4K
Top