1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkRungeKutta2.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 #include "vtkRungeKutta2.h"
16
17 #include "vtkFunctionSet.h"
18 #include "vtkObjectFactory.h"
19
20 vtkStandardNewMacro(vtkRungeKutta2);
21
vtkRungeKutta2()22 vtkRungeKutta2::vtkRungeKutta2()
23 {
24 }
25
~vtkRungeKutta2()26 vtkRungeKutta2::~vtkRungeKutta2()
27 {
28 }
29
30 // Calculate next time step
ComputeNextStep(double * xprev,double * dxprev,double * xnext,double t,double & delT,double & delTActual,double,double,double,double & error)31 int vtkRungeKutta2::ComputeNextStep(double* xprev, double* dxprev, double* xnext,
32 double t, double& delT, double& delTActual,
33 double, double, double, double& error)
34 {
35 int i, numDerivs, numVals;
36
37 delTActual = delT;
38 error = 0.0;
39
40 if (!this->FunctionSet)
41 {
42 vtkErrorMacro("No derivative functions are provided!");
43 return NOT_INITIALIZED;
44 }
45
46 if (!this->Initialized)
47 {
48 vtkErrorMacro("Integrator not initialized!");
49 return NOT_INITIALIZED;
50 }
51
52 numDerivs = this->FunctionSet->GetNumberOfFunctions();
53 numVals = numDerivs + 1;
54 for(i=0; i<numVals-1; i++)
55 {
56 this->Vals[i] = xprev[i];
57 }
58 this->Vals[numVals-1] = t;
59
60 // Obtain the derivatives dx_i at x_i
61 if (dxprev)
62 {
63 for(i=0; i<numDerivs; i++)
64 {
65 this->Derivs[i] = dxprev[i];
66 }
67 }
68 else if ( !this->FunctionSet->FunctionValues(this->Vals, this->Derivs) )
69 {
70 memcpy(xnext, this->Vals, (numVals-1)*sizeof(double));
71 return OUT_OF_DOMAIN;
72 }
73
74 // Half-step
75 for(i=0; i<numVals-1; i++)
76 {
77 this->Vals[i] = xprev[i] + delT/2.0*this->Derivs[i];
78 }
79 this->Vals[numVals-1] = t + delT/2.0;
80
81 // Obtain the derivatives at x_i + dt/2 * dx_i
82 if (!this->FunctionSet->FunctionValues(this->Vals, this->Derivs))
83 {
84 memcpy(xnext, this->Vals, (numVals-1)*sizeof(double));
85 return OUT_OF_DOMAIN;
86 }
87
88 // Calculate x_i using improved values of derivatives
89 for(i=0; i<numDerivs; i++)
90 {
91 xnext[i] = xprev[i] + delT*this->Derivs[i];
92 }
93
94 return 0;
95 }
96
97
98
99
100
101
102