1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
|
/*=========================================================================
Program: Visualization Toolkit
Module: $RCSfile: vtkRungeKutta45.h,v $
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
// .NAME vtkRungeKutta45 - Integrate an initial value problem using 5th
// order Runge-Kutta method with adaptive stepsize control.
// .SECTION Description
// This is a concrete sub-class of vtkInitialValueProblemSolver.
// It uses a 5th order Runge-Kutta method with stepsize control to obtain
// the values of a set of functions at the next time step. The stepsize
// is adjusted by calculating an estimated error using an embedded 4th
// order Runge-Kutta formula:
// Press, W. H. et al., 1992, Numerical Recipes in Fortran, Second
// Edition, Cambridge University Press
// Cash, J.R. and Karp, A.H. 1990, ACM Transactions on Mathematical
// Software, vol 16, pp 201-222
// .SECTION See Also
// vtkInitialValueProblemSolver vtkRungeKutta4 vtkRungeKutta2 vtkFunctionSet
#ifndef __vtkRungeKutta45_h
#define __vtkRungeKutta45_h
#include "vtkInitialValueProblemSolver.h"
class VTK_COMMON_EXPORT vtkRungeKutta45 : public vtkInitialValueProblemSolver
{
public:
vtkTypeRevisionMacro(vtkRungeKutta45,vtkInitialValueProblemSolver);
void PrintSelf(ostream& os, vtkIndent indent);
// Description:
// Construct a vtkRungeKutta45 with no initial FunctionSet.
static vtkRungeKutta45 *New();
// Description:
// Given initial values, xprev , initial time, t and a requested time
// interval, delT calculate values of x at t+delTActual (xnext).
// Possibly delTActual != delT. This may occur
// because this solver supports adaptive stepsize control. It tries
// to change to stepsize such that
// the (estimated) error of the integration is less than maxError.
// The solver will not set the stepsize smaller than minStep or
// larger than maxStep (note that maxStep and minStep should both
// be positive, whereas delT can be negative).
// Also note that delT is an in/out argument. vtkRungeKutta45
// will modify delT to reflect the best (estimated) size for the next
// integration step.
// An estimated value for the error is returned (by reference) in error.
// This is the norm of the error vector if there are more than
// one function to be integrated.
// This method returns an error code representing the nature of
// the failure:
// OutOfDomain = 1,
// NotInitialized = 2,
// UnexpectedValue = 3
virtual int ComputeNextStep(double* xprev, double* xnext, double t,
double& delT, double maxError, double& error)
{
double minStep = delT;
double maxStep = delT;
double delTActual;
return this->ComputeNextStep(xprev, 0, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
}
virtual int ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT,
double maxError, double& error)
{
double minStep = delT;
double maxStep = delT;
double delTActual;
return this->ComputeNextStep(xprev, dxprev, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
}
virtual int ComputeNextStep(double* xprev, double* xnext,
double t, double& delT, double& delTActual,
double minStep, double maxStep,
double maxError, double& error)
{
return this->ComputeNextStep(xprev, 0, xnext, t, delT, delTActual,
minStep, maxStep, maxError, error);
}
virtual int ComputeNextStep(double* xprev, double* dxprev, double* xnext,
double t, double& delT, double& delTActual,
double minStep, double maxStep,
double maxError, double& error);
protected:
vtkRungeKutta45();
~vtkRungeKutta45();
virtual void Initialize();
// Cash-Karp parameters
static double A[5];
static double B[5][5];
static double C[6];
static double DC[6];
double* NextDerivs[6];
int ComputeAStep(double* xprev, double* dxprev, double* xnext, double t,
double& delT, double& error);
private:
vtkRungeKutta45(const vtkRungeKutta45&); // Not implemented.
void operator=(const vtkRungeKutta45&); // Not implemented.
};
#endif
|
