From 7a44617da04bd73f36ae4bc9fc617362d1dfecab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alexander=20D=C3=B6rr?= Date: Tue, 9 Feb 2010 11:09:19 +0000 Subject: [PATCH] fixed bug with event processing --- src/eva2/tools/math/des/RKSolverV2.java | 122 +++++++++++++----------- 1 file changed, 67 insertions(+), 55 deletions(-) diff --git a/src/eva2/tools/math/des/RKSolverV2.java b/src/eva2/tools/math/des/RKSolverV2.java index 963d3826..ae67855e 100644 --- a/src/eva2/tools/math/des/RKSolverV2.java +++ b/src/eva2/tools/math/des/RKSolverV2.java @@ -23,6 +23,7 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable { * */ private static boolean useLinearCalc = true; + /** * * @param args @@ -96,13 +97,14 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable { * @param Ytemp * @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][]; K[0] = Mathematics.svMult(h, EDES.getValue(x, Ytemp)); - K[1] = Mathematics.svMult(h, EDES.getValue(x + h / 2, Mathematics.vvAdd( - Ytemp, Mathematics.svMult(0.5, K[0])))); - K[2] = Mathematics.svMult(h, EDES.getValue(x + h / 2, Mathematics.vvAdd( - Ytemp, Mathematics.svMult(0.5, K[1])))); + K[1] = Mathematics.svMult(h, EDES.getValue(x + h / 2, Mathematics + .vvAdd(Ytemp, Mathematics.svMult(0.5, K[0])))); + K[2] = Mathematics.svMult(h, EDES.getValue(x + h / 2, Mathematics + .vvAdd(Ytemp, Mathematics.svMult(0.5, K[1])))); K[3] = Mathematics.svMult(h, EDES.getValue(x + h, Mathematics.vvAdd( Ytemp, K[2]))); @@ -233,8 +235,8 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable { return solveAtTimePoints(DES, initialValues, timeVector); } - public double[][] solveAtTimePoints(EventDESystem EDES, double[] initialValues, - double[] timePoints) { + public double[][] solveAtTimePoints(EventDESystem EDES, + double[] initialValues, double[] timePoints) { return solveAtTimePoints(EDES, initialValues, timePoints, false); } @@ -257,8 +259,8 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable { /** * */ - public double[][] solveAtTimePointsWithInitialConditions(EventDESystem EDES, - double[][] initConditions, double[] timePoints) { + public double[][] solveAtTimePointsWithInitialConditions( + EventDESystem EDES, double[][] initConditions, double[] timePoints) { int order = EDES.getDESystemDimension(); double[][] result = new double[timePoints.length][order]; result[0] = initConditions[0]; @@ -302,8 +304,8 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable { * @param timeEnd * @return */ - public double[][] solveByStepSize(EventDESystem EDES, double[] initialValues, - double timeBegin, double timeEnd) { + public double[][] solveByStepSize(EventDESystem EDES, + double[] initialValues, double timeBegin, double timeEnd) { 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); } - public double[] processEvents(double time, double[] Ytemp, EventDESystem EDES) { + public double[] processEvents(double time, double[] Ytemp, + EventDESystem EDES) { double[] res = new double[Ytemp.length]; - //Arrays.fill(res, 0); - double[] delays = EDES.processEvents(time, Ytemp, res); - for (int j = 0; j < delays.length; j++) - if (delays[j] != Double.NaN) { - //System.out.printf("Ytemp[j]_old = %s\tYtemp[j]_new = %s\n", Ytemp[j], res[j]); - Ytemp[j] = res[j]; - } + // Arrays.fill(res, 0); + double[] delays = EDES.processEvents(time, Ytemp, res); + System.out.println("delay " +Arrays.toString(delays)); + for (int j = 0; j < delays.length; j++) + if (!Double.isNaN(delays[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; - + } - - public void processEventsVoid(double time, double[] Ytemp, EventDESystem EDES) { + public static boolean containsNaN(double[] arr){ + 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]; - //Arrays.fill(res, 0); - double[] delays = EDES.processEvents(time, Ytemp, res); - for (int j = 0; j < delays.length; j++) - if (delays[j] != Double.NaN) { - //System.out.printf("Ytemp[j]_old = %s\tYtemp[j]_new = %s\n", Ytemp[j], res[j]); - Ytemp[j] = res[j]; + // Arrays.fill(res, 0); + double[] delays = EDES.processEvents(time, Ytemp, res); + for (int j = 0; j < delays.length; j++) + if (!Double.isNaN(delays[j])) { + // System.out.printf("Ytemp[j]_old = %s\tYtemp[j]_new = %s\n", + // Ytemp[j], res[j]); + Ytemp[j] = res[j]; } - + } + /** * When set to TRUE, includeTimes will make the * 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 */ - private double[][] solveAtTimePoints(EventDESystem EDES, double[] initialValues, - double[] timePoints, boolean includeTimes) { + private double[][] solveAtTimePoints(EventDESystem EDES, + double[] initialValues, double[] timePoints, boolean includeTimes) { // sorted timepoints!!!!!!!!!!!!!!!!!!!!! int order = EDES.getDESystemDimension(); double result[][], x = timePoints[0]; @@ -389,28 +406,21 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable { result[i - 1].length - 1); else Ytemp = result[i - 1].clone(); - - //process events - /*double[] res = new double[Ytemp.length]; - Arrays.fill(res, 0); - 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); - + + // process events + //double[] YtempClone = processEvents(x, Ytemp.clone(), EDES); + // Ytemp = processEventsVoid(x, Ytemp,EDES); + for (int j = 0; j < inbetweensteps; j++) { + if (useLinearCalc) rkTerm2(EDES, h, x, Ytemp, change); else change = rkTerm(EDES, h, x, Ytemp); // System.out.println("aft change 0 " + change[0]); - + Mathematics.vvAdd(Ytemp, change, Ytemp); + if (this.nonnegative) { for (int k = 0; k < Ytemp.length; k++) { @@ -420,12 +430,14 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable { } x += h; - - //double[] YtempClone2 = processEvents(x, Ytemp.clone(),EDES); - processEventsVoid(x, Ytemp,EDES); - //Mathematics.vvAdd(Ytemp, YtempClone2, Ytemp); - + + // process events + double[] YtempClone2 = processEvents(x, Ytemp.clone(), EDES); + Mathematics.vvAdd(Ytemp, YtempClone2, Ytemp); + //processEventsVoid(x, Ytemp,EDES); + } + // ohne wirkung //Mathematics.vvAdd(Ytemp, YtempClone, Ytemp); h = timePoints[i] - x; @@ -436,8 +448,7 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable { change = rkTerm(EDES, h, x, Ytemp); Mathematics.vvAdd(Ytemp, change, Ytemp); - - + if (this.nonnegative) { for (int k = 0; k < Ytemp.length; k++) { if (Ytemp[k] < 0) @@ -464,8 +475,9 @@ public class RKSolverV2 extends AbstractDESSolver implements Serializable { * @param timeEnd * @return */ - private double[][] solveByStepSize(EventDESystem EDES, double[] initialValues, - double timeBegin, double timeEnd, boolean time) { + private double[][] solveByStepSize(EventDESystem EDES, + double[] initialValues, double timeBegin, double timeEnd, + boolean time) { int numsteps = (int) Math.round(((timeEnd - timeBegin) / stepSize) + 1); unstableFlag = false; // System.out.println(numsteps);