halfdan/eva2
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

fixed bug with event processing

Browse Source
This commit is contained in:
Alexander Dörr 2010-02-09 11:09:19 +00:00
parent 41b4f9a322
commit 7a44617da0
1 changed files with 67 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

122
src/eva2/tools/math/des/RKSolverV2.java
View File

@ -23,6 +23,7 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable {
* *
*/ */
private static boolean useLinearCalc = true; private static boolean useLinearCalc = true;
/** /**
* *
* @param args * @param args
@ -96,13 +97,14 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable {
* @param Ytemp * @param Ytemp
* @return * @return
*/ */
public double[] rkTerm(EventDESystem EDES, double h, double x, double[] Ytemp) { public double[] rkTerm(EventDESystem EDES, double h, double x,
double[] Ytemp) {
double[][] K = new double[4][]; double[][] K = new double[4][];
K[0] = Mathematics.svMult(h, EDES.getValue(x, Ytemp)); K[0] = Mathematics.svMult(h, EDES.getValue(x, Ytemp));
K[1] = Mathematics.svMult(h, EDES.getValue(x + h / 2, Mathematics.vvAdd( K[1] = Mathematics.svMult(h, EDES.getValue(x + h / 2, Mathematics
Ytemp, Mathematics.svMult(0.5, K[0])))); .vvAdd(Ytemp, Mathematics.svMult(0.5, K[0]))));
K[2] = Mathematics.svMult(h, EDES.getValue(x + h / 2, Mathematics.vvAdd( K[2] = Mathematics.svMult(h, EDES.getValue(x + h / 2, Mathematics
Ytemp, Mathematics.svMult(0.5, K[1])))); .vvAdd(Ytemp, Mathematics.svMult(0.5, K[1]))));
K[3] = Mathematics.svMult(h, EDES.getValue(x + h, Mathematics.vvAdd( K[3] = Mathematics.svMult(h, EDES.getValue(x + h, Mathematics.vvAdd(
Ytemp, K[2]))); Ytemp, K[2])));
@ -233,8 +235,8 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable {
return solveAtTimePoints(DES, initialValues, timeVector); return solveAtTimePoints(DES, initialValues, timeVector);
} }
public double[][] solveAtTimePoints(EventDESystem EDES, double[] initialValues, public double[][] solveAtTimePoints(EventDESystem EDES,
double[] timePoints) { double[] initialValues, double[] timePoints) {
return solveAtTimePoints(EDES, initialValues, timePoints, false); return solveAtTimePoints(EDES, initialValues, timePoints, false);
} }
@ -257,8 +259,8 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable {
/** /**
* *
*/ */
public double[][] solveAtTimePointsWithInitialConditions(EventDESystem EDES, public double[][] solveAtTimePointsWithInitialConditions(
double[][] initConditions, double[] timePoints) { EventDESystem EDES, double[][] initConditions, double[] timePoints) {
int order = EDES.getDESystemDimension(); int order = EDES.getDESystemDimension();
double[][] result = new double[timePoints.length][order]; double[][] result = new double[timePoints.length][order];
result[0] = initConditions[0]; result[0] = initConditions[0];
@ -302,8 +304,8 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable {
* @param timeEnd * @param timeEnd
* @return * @return
*/ */
public double[][] solveByStepSize(EventDESystem EDES, double[] initialValues, public double[][] solveByStepSize(EventDESystem EDES,
double timeBegin, double timeEnd) { double[] initialValues, double timeBegin, double timeEnd) {
return solveByStepSize(EDES, initialValues, timeBegin, timeEnd, false); return solveByStepSize(EDES, initialValues, timeBegin, timeEnd, false);
} }
@ -319,31 +321,46 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable {
return solveByStepSize(EDES, initialValues, timeBegin, timeEnd, true); return solveByStepSize(EDES, initialValues, timeBegin, timeEnd, true);
} }
public double[] processEvents(double time, double[] Ytemp, EventDESystem EDES) { public double[] processEvents(double time, double[] Ytemp,
EventDESystem EDES) {
double[] res = new double[Ytemp.length]; double[] res = new double[Ytemp.length];
//Arrays.fill(res, 0); // Arrays.fill(res, 0);
double[] delays = EDES.processEvents(time, Ytemp, res); double[] delays = EDES.processEvents(time, Ytemp, res);
for (int j = 0; j < delays.length; j++) System.out.println("delay " +Arrays.toString(delays));
if (delays[j] != Double.NaN) { for (int j = 0; j < delays.length; j++)
//System.out.printf("Ytemp[j]_old = %s\tYtemp[j]_new = %s\n", Ytemp[j], res[j]); if (!Double.isNaN(delays[j])) {
Ytemp[j] = res[j]; System.out.printf("time %s :Ytemp[%s]_old = %s\tYtemp[%s]_new = %s\n",
} time,j, Ytemp[j], j, (res[j] - Ytemp[j]));
Ytemp[j] = res[j] - Ytemp[j];
} else Ytemp[j] = 0d;
return Ytemp; return Ytemp;
} }
public static boolean containsNaN(double[] arr){
public void processEventsVoid(double time, double[] Ytemp, EventDESystem EDES) { for (int i = 0; i < arr.length; i++) {
if (Double.isNaN(arr[i]))
return true;
}
return false;
}
public void processEventsVoid(double time, double[] Ytemp,
EventDESystem EDES) {
double[] res = new double[Ytemp.length]; double[] res = new double[Ytemp.length];
//Arrays.fill(res, 0); // Arrays.fill(res, 0);
double[] delays = EDES.processEvents(time, Ytemp, res); double[] delays = EDES.processEvents(time, Ytemp, res);
for (int j = 0; j < delays.length; j++) for (int j = 0; j < delays.length; j++)
if (delays[j] != Double.NaN) { if (!Double.isNaN(delays[j])) {
//System.out.printf("Ytemp[j]_old = %s\tYtemp[j]_new = %s\n", Ytemp[j], res[j]); // System.out.printf("Ytemp[j]_old = %s\tYtemp[j]_new = %s\n",
Ytemp[j] = res[j]; // Ytemp[j], res[j]);
Ytemp[j] = res[j];
} }
} }
/** /**
* When set to <code>TRUE</code>, <code>includeTimes</code> will make the * When set to <code>TRUE</code>, <code>includeTimes</code> will make the
* solver to return a matrix with the first column containing the times. By * solver to return a matrix with the first column containing the times. By
@ -351,8 +368,8 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable {
* *
* @param includeTimes * @param includeTimes
*/ */
private double[][] solveAtTimePoints(EventDESystem EDES, double[] initialValues, private double[][] solveAtTimePoints(EventDESystem EDES,
double[] timePoints, boolean includeTimes) { double[] initialValues, double[] timePoints, boolean includeTimes) {
// sorted timepoints!!!!!!!!!!!!!!!!!!!!! // sorted timepoints!!!!!!!!!!!!!!!!!!!!!
int order = EDES.getDESystemDimension(); int order = EDES.getDESystemDimension();
double result[][], x = timePoints[0]; double result[][], x = timePoints[0];
@ -389,28 +406,21 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable {
result[i - 1].length - 1); result[i - 1].length - 1);
else else
Ytemp = result[i - 1].clone(); Ytemp = result[i - 1].clone();
//process events // process events
/*double[] res = new double[Ytemp.length]; //double[] YtempClone = processEvents(x, Ytemp.clone(), EDES);
Arrays.fill(res, 0); // Ytemp = processEventsVoid(x, Ytemp,EDES);
double[] delays = EDES.processEvents(x, Ytemp, res);
for (int j = 0; j < delays.length; j++)
if (delays[j] != Double.NaN) {
System.out.println(res[j]);
Ytemp[j] = res[j];
}
*/
//double[] YtempClone = processEvents(x, Ytemp.clone(),EDES);
processEventsVoid(x, Ytemp,EDES);
for (int j = 0; j < inbetweensteps; j++) { for (int j = 0; j < inbetweensteps; j++) {
if (useLinearCalc) if (useLinearCalc)
rkTerm2(EDES, h, x, Ytemp, change); rkTerm2(EDES, h, x, Ytemp, change);
else else
change = rkTerm(EDES, h, x, Ytemp); change = rkTerm(EDES, h, x, Ytemp);
// System.out.println("aft change 0 " + change[0]); // System.out.println("aft change 0 " + change[0]);
Mathematics.vvAdd(Ytemp, change, Ytemp); Mathematics.vvAdd(Ytemp, change, Ytemp);
if (this.nonnegative) { if (this.nonnegative) {
for (int k = 0; k < Ytemp.length; k++) { for (int k = 0; k < Ytemp.length; k++) {
@ -420,12 +430,14 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable {
} }
x += h; x += h;
//double[] YtempClone2 = processEvents(x, Ytemp.clone(),EDES); // process events
processEventsVoid(x, Ytemp,EDES); double[] YtempClone2 = processEvents(x, Ytemp.clone(), EDES);
//Mathematics.vvAdd(Ytemp, YtempClone2, Ytemp); Mathematics.vvAdd(Ytemp, YtempClone2, Ytemp);
//processEventsVoid(x, Ytemp,EDES);
} }
// ohne wirkung
//Mathematics.vvAdd(Ytemp, YtempClone, Ytemp); //Mathematics.vvAdd(Ytemp, YtempClone, Ytemp);
h = timePoints[i] - x; h = timePoints[i] - x;
@ -436,8 +448,7 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable {
change = rkTerm(EDES, h, x, Ytemp); change = rkTerm(EDES, h, x, Ytemp);
Mathematics.vvAdd(Ytemp, change, Ytemp); Mathematics.vvAdd(Ytemp, change, Ytemp);
if (this.nonnegative) { if (this.nonnegative) {
for (int k = 0; k < Ytemp.length; k++) { for (int k = 0; k < Ytemp.length; k++) {
if (Ytemp[k] < 0) if (Ytemp[k] < 0)
@ -464,8 +475,9 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable {
* @param timeEnd * @param timeEnd
* @return * @return
*/ */
private double[][] solveByStepSize(EventDESystem EDES, double[] initialValues, private double[][] solveByStepSize(EventDESystem EDES,
double timeBegin, double timeEnd, boolean time) { double[] initialValues, double timeBegin, double timeEnd,
boolean time) {
int numsteps = (int) Math.round(((timeEnd - timeBegin) / stepSize) + 1); int numsteps = (int) Math.round(((timeEnd - timeBegin) / stepSize) + 1);
unstableFlag = false; unstableFlag = false;
// System.out.println(numsteps); // System.out.println(numsteps);