diff --git a/tests/TestHelperFunctions.m b/tests/TestHelperFunctions.m index 93c17dcd..30da0c7b 100644 --- a/tests/TestHelperFunctions.m +++ b/tests/TestHelperFunctions.m @@ -26,7 +26,7 @@ function helperFunctionInElseBlock(this) datahandle = prepareDatahandleForIntegration( ... 'funinelse_RHS', ... - 'integrator', func2str(integrator), ... + 'integrator', integrator, ... 'options', odeoptions); tEnd = 20; tspan = [0 tEnd]; @@ -44,7 +44,7 @@ function manyFunctions(this) datahandle = prepareDatahandleForIntegration( ... 'manyFunctions', ... - 'integrator', func2str(integrator), ... + 'integrator', integrator, ... 'options', odeoptions); tEnd = 30; tspan = [0 tEnd]; diff --git a/tests/TestIgnore.m b/tests/TestIgnore.m index 599506f3..0d8a5ee4 100644 --- a/tests/TestIgnore.m +++ b/tests/TestIgnore.m @@ -21,7 +21,7 @@ function testIgnoreSimple(testCase) p = 0; datahandle = prepareDatahandleForIntegration('simpleIgnoreRHS', ... - 'solver', func2str(solver), ... + 'solver', solver, ... 'options', options); data = datahandle.getData(); % ensure only one of the two ifs got turned into a ctrlif @@ -42,7 +42,7 @@ function testIgnoreNested(testCase) p = 0; datahandle = prepareDatahandleForIntegration('nestedIgnoreRHS', ... - 'solver', func2str(solver), ... + 'solver', solver, ... 'options', options); data = datahandle.getData(); % ensure only one of the two ifs got turned into a ctrlif @@ -63,7 +63,7 @@ function testIgnoredIfInNonIgnoredIf(testCase) p = 0; datahandle = prepareDatahandleForIntegration('ignoredIfInNonIgnoredIfRHS', ... - 'solver', func2str(solver), ... + 'solver', solver, ... 'options', options); data = datahandle.getData(); % ensure only the outer if got turned into a ctrlif @@ -84,7 +84,7 @@ function testIgnoredIfInHelper(testCase) p = 0; datahandle = prepareDatahandleForIntegration('ignoreInHelperRHS', ... - 'solver', func2str(solver), ... + 'solver', solver, ... 'options', options); data = datahandle.getData(); % ensure only one of the two ifs got turned into a ctrlif @@ -108,7 +108,7 @@ function testHelperInIgnore(testCase) p = 0; datahandle = prepareDatahandleForIntegration('helperInIgnoreRHS', ... - 'solver', func2str(solver), ... + 'solver', solver, ... 'options', options); data = datahandle.getData(); @@ -131,7 +131,7 @@ function testIgnoreAndNotIgnore(testCase) p = 0; datahandle = prepareDatahandleForIntegration('ignoreAndNotIgnoreRHS', ... - 'solver', func2str(solver), ... + 'solver', solver, ... 'options', options); data = datahandle.getData(); diff --git a/tests/TestSensitivities.m b/tests/TestSensitivities.m index 0d06e225..f7a73409 100644 --- a/tests/TestSensitivities.m +++ b/tests/TestSensitivities.m @@ -11,8 +11,8 @@ function testCanonexDirectional(testCase, method) [integrator, options, t0, tEnd, p, x0] = getOdeDataForCanonex(testCase); datahandle = prepareDatahandleForIntegration( ... - 'canonicalExampleRHS', ... - 'solver', func2str(integrator), ... + 'rhsCanonicalExample', ... + 'solver', integrator, ... 'options', options); sol = solveODE(datahandle, [t0 tEnd], x0, p); @@ -34,8 +34,8 @@ function testCanonexVde(testCase) % Test the sensitivities generated with the VDE method on the canonical example. [integrator, options, t0, tEnd, p, x0] = getOdeDataForCanonex(testCase); datahandle = prepareDatahandleForIntegration( ... - 'canonicalExampleRHS', ... - 'solver', func2str(integrator), ... + 'rhsCanonicalExample', ... + 'solver', integrator, ... 'options', options); sol = solveODE(datahandle, [t0 tEnd], x0, p); @@ -80,7 +80,7 @@ function testCanonexEndPiecewise(testCase) % Test the sensitivities generated with the END_piecewise method on the canonical example. [integrator, options, t0, tEnd, p, x0] = getOdeDataForCanonex(testCase); datahandle = prepareDatahandleForIntegration( ... - 'canonicalExampleRHS', ... + 'rhsCanonicalExample', ... 'solver', func2str(integrator), ... 'options', options); sol = solveODE(datahandle, [t0 tEnd], x0, p); @@ -129,7 +129,7 @@ function testCanonexEndFull(testCase) % instead of ts+eps(ts). [integrator, options, t0, tEnd, p, x0] = getOdeDataForCanonex(testCase); datahandle = prepareDatahandleForIntegration( ... - 'canonicalExampleRHS', ... + 'rhsCanonicalExample', ... 'solver', func2str(integrator), ... 'options', options); sol = solveODE(datahandle, [t0 tEnd], x0, p); diff --git a/tests/TestSolverCompatibility.m b/tests/TestSolverCompatibility.m index f00cb997..27caa6ce 100644 --- a/tests/TestSolverCompatibility.m +++ b/tests/TestSolverCompatibility.m @@ -21,7 +21,7 @@ function setTimer(testCase) canonex_tspan = [0 20]; canonex_x0 = [1;0]; canonex_p = 5.437; - canonex_rhsFunction = 'canonicalExampleRHS'; + canonex_rhsFunction = 'rhsCanonicalExample'; canonex_odeoptions = odeset('AbsTol', 1e-12, 'RelTol', 1e-10); canonex_xEnd = [49.255067; 1.348824]; canonex_SWPs = [11.000275475 11.27004032]; @@ -29,19 +29,19 @@ function setTimer(testCase) subway_tspan = [0 65]; subway_x0 = [0; 0; 0]; % subway_p initialized later because we need to initIFDIFF() before - subway_rhsFunction = 'newYorkCitySubwayModelRhs'; + subway_rhsFunction = 'rhsNewYorkCitySubwayModel.m'; subway_odeoptions = odeset( 'AbsTol', 1e-20, 'RelTol', 1e-10); subway_xEnd = [2112.07361577; 0.00124794; 4124.77885608]; subway_SWPs = [0.63166061, 2.43955402, 3.64338000, 5.60010643, 12.60705000, 45.78275000, 57.16005000]; % DAE example: input parameters and expected results - dae_tspan = [0 5]; + dae_tspan = [0 3]; + dae_p = -0.3 dae_x0 = [1; -1]; - dae_p = -0.2 - dae_rhsFunction = 'daeExampleRHS'; - dae_odeoptions = odeset('Mass', [1 0; 0 0], 'MassSingular', 'yes', 'AbsTol', 1e-9,'RelTol', 1e-6); - dae_xEnd = [0.199999951200446; -0.199999951200446] - dae_SWPs = [1.609435443361949]; + dae_rhsFunction = 'rhsDaeExampleWorkshop'; + dae_odeoptions = odeset('Mass', [1 0; 0 0], 'MassSingular', 'yes', 'AbsTol', 1e-9,'RelTol', 1e-5); + dae_xEnd = [0.427392826365841; -0.299999535808094] + dae_SWPs = [1.068023313185268]; end properties (Access = private) @@ -185,8 +185,8 @@ function testOde23tDAE(testCase) % lower tolerance since ode23t is an order 2 method % (uses trapezoidal rule) - testCase.verifyEqual(sol.y(:, end), testCase.expected_xEnd, "RelTol", 1e-5); - testCase.verifyEqual(sol.switches, testCase.expected_SWPs, "RelTol", 1e-4); + testCase.verifyEqual(sol.y(:, end), testCase.expected_xEnd, "RelTol", 1e-3); + testCase.verifyEqual(sol.switches, testCase.expected_SWPs, "RelTol", 1e-3); end function testOde23tbCanonex(testCase) setCanonexParameters(testCase); diff --git a/tests/TestStateJumps.m b/tests/TestStateJumps.m index 2da8e876..8f094a68 100644 --- a/tests/TestStateJumps.m +++ b/tests/TestStateJumps.m @@ -28,7 +28,7 @@ function testIdenticalIfs(testCase) options = odeset('AbsTol', 1e-8, 'RelTol', 1e-6); datahandle = prepareDatahandleForIntegration( ... 'identicalIfRHS', ... - 'solver', func2str(integrator), ... + 'solver', integrator, ... 'options', options); t0 = 0; tEnd = 3; @@ -57,7 +57,7 @@ function testTwoJumpsOppositeDirections(testCase) options = odeset('AbsTol', 1e-8, 'RelTol', 1e-6); datahandle = prepareDatahandleForIntegration( ... 'twoJumpsGoodRHS', ... - 'solver', func2str(integrator), ... + 'solver', integrator, ... 'options', options); t0 = 0; tEnd = 4.1; % switches at tEnd exactly unfortunately aren't handled too gracefully @@ -88,7 +88,7 @@ function testTwoJumpsSameDirection(testCase) datahandle = prepareDatahandleForIntegration( ... 'twoJumpsBadRHS', ... - 'solver', func2str(integrator), ... + 'solver', integrator, ... 'options', options); testCase.verifyError(@() solveODE(datahandle, [t0 tF], x0, p), ''); @@ -100,7 +100,7 @@ function testJumpChangesModel(testCase) options = odeset('AbsTol', 1e-8, 'RelTol', 1e-6); datahandle = prepareDatahandleForIntegration( ... 'jumpChangesModelRHS', ... - 'solver', func2str(integrator), ... + 'solver', integrator, ... 'options', options); t0 = 0; tEnd = 4; @@ -118,7 +118,7 @@ function testJumpInHelper(testCase) options = odeset('AbsTol', 1e-8, 'RelTol', 1e-6); datahandle = prepareDatahandleForIntegration( ... 'jumpInHelperRHS', ... - 'solver', func2str(integrator), ... + 'solver', integrator, ... 'options', options); t0 = 0; tEnd = 25; @@ -146,7 +146,7 @@ function testJumpInHelperWithJumpsDisabled(testCase) options = odeset('AbsTol', 1e-8, 'RelTol', 1e-6); datahandle = prepareDatahandleForIntegration( ... 'jumpInHelperRHS', ... - 'solver', func2str(integrator), ... + 'solver', integrator, ... 'options', options); t0 = 0; tEnd = 25; @@ -171,8 +171,8 @@ function testBounceball(testCase) integrator = TestStateJumps.defaultIntegrator; options = odeset('AbsTol', 1e-8, 'RelTol', 1e-6); datahandle = prepareDatahandleForIntegration( ... - 'bounceballRHS', ... - 'solver', func2str(integrator), ... + 'rhsBounceball', ... + 'solver', integrator, ... 'options', options); g = 9.807; gamma = 0.9; @@ -230,7 +230,7 @@ function testTwoUpdates(testCase) options = odeset('AbsTol', 1e-8, 'RelTol', 1e-6); datahandle = prepareDatahandleForIntegration( ... 'twoUpdatesRHS', ... - 'solver', func2str(integrator), ... + 'solver', integrator, ... 'options', options); sol = solveODE(datahandle, [0 7], 0, 0); testCase.verifyEqual(length(sol.switches), 2); @@ -241,13 +241,13 @@ function testTwoUpdates(testCase) end function testSensitivitiesSimpleVDE(testCase) - % Test sensitivity computation across jumps using the simple, one-dimensional jumpSensitivityRHS. + % Test sensitivity computation across jumps using the simple, one-dimensional rhsJumpSensitivity. integrator = TestStateJumps.defaultIntegrator; options = odeset('AbsTol', 1e-8, 'RelTol', 1e-6); [t0, tEnd, p, x0] = jumpSensitivityInitdata(); datahandle = prepareDatahandleForIntegration( ... - 'jumpSensitivityRHS', ... - 'solver', func2str(integrator), ... + 'rhsJumpSensitivity', ... + 'solver', integrator, ... 'options', options); sol = solveODE(datahandle, [t0 tEnd], x0, p); testCase.verifyEqual(length(sol.switches), 1); @@ -277,13 +277,13 @@ function testSensitivitiesSimpleVDE(testCase) testCase.verifyEqual(Gp{4}, Gy2(tEnd) * (Uy1 * Gp1(t1Minus) + Up1) + Gp2(tEnd), 'RelTol', rtol2); end function testSensitivitiesSimpleEND_piecewise(testCase) - % Test sensitivity computation across jumps using the simple, one-dimensional jumpSensitivityRHS. + % Test sensitivity computation across jumps using the simple, one-dimensional rhsJumpSensitivity. integrator = TestStateJumps.defaultIntegrator; options = odeset('AbsTol', 1e-8, 'RelTol', 1e-6); [t0, tEnd, p, x0] = jumpSensitivityInitdata(); datahandle = prepareDatahandleForIntegration( ... - 'jumpSensitivityRHS', ... - 'solver', func2str(integrator), ... + 'rhsJumpSensitivity', ... + 'solver', integrator, ... 'options', options); sol = solveODE(datahandle, [t0 tEnd], x0, p); testCase.verifyEqual(length(sol.switches), 1); @@ -314,13 +314,13 @@ function testSensitivitiesSimpleEND_piecewise(testCase) testCase.verifyEqual(Gp{4}, Gy2(tEnd) * (Uy1 * Gp1(t1Minus) + Up1) + Gp2(tEnd), 'RelTol', rtol2); end function testSensitivitiesSimpleEND_full(testCase) - % Test sensitivity computation across jumps using the simple, one-dimensional jumpSensitivityRHS. + % Test sensitivity computation across jumps using the simple, one-dimensional rhsJumpSensitivity. integrator = TestStateJumps.defaultIntegrator; options = odeset('AbsTol', 1e-8, 'RelTol', 1e-6); [t0, tEnd, p, x0] = jumpSensitivityInitdata(); datahandle = prepareDatahandleForIntegration( ... - 'jumpSensitivityRHS', ... - 'solver', func2str(integrator), ... + 'rhsJumpSensitivity', ... + 'solver', integrator, ... 'options', options); sol = solveODE(datahandle, [t0 tEnd], x0, p); testCase.verifyEqual(length(sol.switches), 1); @@ -367,7 +367,7 @@ function testHelperInJump(testCase) datahandle = prepareDatahandleForIntegration( ... 'helperInJumpRHS', ... - 'solver', func2str(integrator), ... + 'solver', integrator, ... 'options', options); sol = solveODE(datahandle, [t0 tF], x0, p); @@ -386,7 +386,7 @@ function testHelperInJumpInHelper(testCase) options = odeset('AbsTol', 1e-8, 'RelTol', 1e-6); datahandle = prepareDatahandleForIntegration( ... 'helperInJumpInHelperRHS', ... - 'solver', func2str(integrator), ... + 'solver', integrator, ... 'options', options); t0 = 0; tEnd = 3; diff --git a/toolbox/examples/bounceball/bounceball_rhs.m b/toolbox/examples/bounceball/bounceball_rhs.m deleted file mode 100644 index 72a94a62..00000000 --- a/toolbox/examples/bounceball/bounceball_rhs.m +++ /dev/null @@ -1,28 +0,0 @@ -function dx = bounceball_rhs(t, x, p) - % t: time - % x: state - % p: parameter vector - % p(1) is a bounce-parameter (currently not in use!) - % p(2) is object's mass - % p(3) is gravitational force - % - % This example implements a simple model of a ball bouncing in - % one dimension, i.e., horizontally. - % To avoid state-jumps, we allow the position x(1) of the ball - % to be negative. After integration, one should take the absolute - % value of x(1). - - dx = zeros(2,1); - m = p(2); - - % RHS logic - if x(1)>=0 - g = exp(2*t)*p(3); - else - g = -exp(2*t)*p(3); - end - - dx(1) = x(2); - dx(2) = -g/m; - -end \ No newline at end of file diff --git a/toolbox/examples/bounceball/bounceball_test.m b/toolbox/examples/bounceball/bounceball_test.m deleted file mode 100644 index 77849414..00000000 --- a/toolbox/examples/bounceball/bounceball_test.m +++ /dev/null @@ -1,36 +0,0 @@ -integrator = @ode45; -t0 = 0; -tf = 3; -timeinterval = [t0,tf]; -initstates = [0, 1]; -p = [0.9, 1, 9.81]; - -fprintf('Preprocessing...\n '); -odeoptions = odeset( 'AbsTol', 1e-14, 'RelTol', 1e-12); -filename = 'bounceball_rhs'; -tic -handlerhs_test = prepareDatahandleForIntegration(filename, 'integrator', func2str(integrator), 'options', odeoptions); -toc - -fprintf('Integration with ifdiff/%s...\n ', func2str(integrator)); -tic -sol_rhs_test = solveODE(handlerhs_test, timeinterval, initstates, p); -toc - -% Switches -sw_pts = cell2mat(handlerhs_test.getData().SWP_detection.switchingpoints); -sw_pts_y = deval(sol_rhs_test, sw_pts); -n_sw_pts = max(size(handlerhs_test.getData().SWP_detection.switchingpoints)); -fprintf('Total number of switches during integration: %d \n', n_sw_pts); - -% Visualize solution -clf('reset'); -axis([0 20 0 1.1]); -lineWidth = 0.8; - -T_gridsize = 0.0001; -T = t0:T_gridsize:tf; -Y = deval(sol_rhs_test,T); -plot(T, abs(Y(1,:)), sw_pts, abs(sw_pts_y(1,:)), 'o', 'lineWidth', lineWidth); -legend('x_1(t)','Switches', 'Location', 'North'); -sgtitle("Zeno's phenomenon in model of bouncing ball"); \ No newline at end of file diff --git a/toolbox/examples/canonicalExample/canonicalExampleRHS_extended.m b/toolbox/examples/canonicalExample/canonicalExampleRHS_extended.m deleted file mode 100644 index 9b76bf59..00000000 --- a/toolbox/examples/canonicalExample/canonicalExampleRHS_extended.m +++ /dev/null @@ -1,32 +0,0 @@ -function dx = canonicalExampleRHS_extended(t,x,p) -dx = zeros(2,1); -%dx(1) = 0.01 * t.^2 + x(2).^2; -dx(1) = 0.01*t + x(2); -dx(2) = 0; - -% a = 1:10:41; -% dx(2) = aNiceFunction(x,t,p, a, dx(2)); - -% dx(2) = checkFeasibility(... -% x, t, p, ... % standardinput (in general for switchingpoint detection not necessarly required) -% 2, 3, ... % intervalboarders -% 1, dx(2), 1); % - - -a = 1:10:101; - -dx(2) = aNiceFunction(x,t,p,a, dx(2)); - -a = a + 5000; -dx(2) = aNiceFunction(x,t,p,a, dx(2)); - -a = a + 5000; -dx(2) = aNiceFunction(x,t,p,a, dx(2)); - -a = a + 5000; -dx(2) = aNiceFunction(x,t,p,a, dx(2)); - - - -end - diff --git a/toolbox/examples/canonicalExample/canonicalExampleRHS_extended_test.m b/toolbox/examples/canonicalExample/canonicalExampleRHS_extended_test.m deleted file mode 100644 index 5eae7939..00000000 --- a/toolbox/examples/canonicalExample/canonicalExampleRHS_extended_test.m +++ /dev/null @@ -1,40 +0,0 @@ -integrator = @ode45; -timeinterval = [0, 100]; -initstates = [1 0]; -%p = 5.437; -p = 0.5; -%% - -options = odeset('AbsTol', 1e-14, 'RelTol', 1e-12); - -filename = 'canonicalExampleRHS_extended'; -clc -tic -hdlrhs_test = prepareDatahandleForIntegration(filename, 'integrator', func2str(integrator), 'options', options); -toc - - -%% - -clc - -tic -[sol_rhs_test, datanyc] = solveODE(hdlrhs_test, timeinterval, initstates, p); -toc - -%% -tic -ode = ode45(@(t,x) canonicalExampleRHS_extended(t,x,p), timeinterval, initstates, options); -toc - - -%% -hold on -plot(sol_rhs_test.x, sol_rhs_test.y(1,:), '*r', sol_rhs_test.x, sol_rhs_test.y(2,:), '*r') -plot(ode.x, ode.y(1,:), 'k', ode.x, ode.y(2,:), 'k') -hold off - - - - - diff --git a/toolbox/examples/canonicalExample/canonicalExampleRHS_reformulated.m b/toolbox/examples/canonicalExample/canonicalExampleRHS_reformulated.m deleted file mode 100644 index df405751..00000000 --- a/toolbox/examples/canonicalExample/canonicalExampleRHS_reformulated.m +++ /dev/null @@ -1,19 +0,0 @@ -function dx = canonicalExampleRHS_reformulated(t,x,p) - -dx = zeros(2,1); -dx(1) = 0.01 * t.^2 + x(2).^3; - -if x(1) > p(1) + 0.5 - dx(2) = 0; -else - if x(1) > p(1) - dx(2) = 5; - else - dx(2) = 0; - end -end - -end - - - diff --git a/toolbox/examples/canonicalExample/canonicalExample_naive.m b/toolbox/examples/canonicalExample/naiveCanonicalExample.m similarity index 86% rename from toolbox/examples/canonicalExample/canonicalExample_naive.m rename to toolbox/examples/canonicalExample/naiveCanonicalExample.m index 730dce7d..c6cda68e 100644 --- a/toolbox/examples/canonicalExample/canonicalExample_naive.m +++ b/toolbox/examples/canonicalExample/naiveCanonicalExample.m @@ -13,11 +13,11 @@ % integrate fprintf('Integrating ...') tic -sol = integrator(@(t,x) canonicalExampleRHS(t,x,p), tspan, y0, options); +sol = integrator(@(t,x) rhsCanonicalExample(t,x,p), tspan, y0, options); toc fprintf('Integrating with limited step size...') tic -soltrue = integrator(@(t,x) canonicalExampleRHS(t,x,p), tspan, y0, optionsExact); +soltrue = integrator(@(t,x) rhsCanonicalExample(t,x,p), tspan, y0, optionsExact); toc % evaluate solution diff --git a/toolbox/examples/canonicalExample/canonicalExampleRHS.m b/toolbox/examples/canonicalExample/rhsCanonicalExample.m similarity index 81% rename from toolbox/examples/canonicalExample/canonicalExampleRHS.m rename to toolbox/examples/canonicalExample/rhsCanonicalExample.m index 28f0344d..49996635 100644 --- a/toolbox/examples/canonicalExample/canonicalExampleRHS.m +++ b/toolbox/examples/canonicalExample/rhsCanonicalExample.m @@ -1,4 +1,4 @@ -function dx = canonicalExampleRHS(t,x,p) +function dx = rhsCanonicalExample(t,x,p) dx = zeros(2,1); dx(1) = 0.01 * t.^2 + x(2).^3; diff --git a/toolbox/examples/canonicalExample/analyticalSolution.m b/toolbox/examples/canonicalExample/runAnalyticalSolution.m similarity index 91% rename from toolbox/examples/canonicalExample/analyticalSolution.m rename to toolbox/examples/canonicalExample/runAnalyticalSolution.m index b3461fb9..cebac8a5 100644 --- a/toolbox/examples/canonicalExample/analyticalSolution.m +++ b/toolbox/examples/canonicalExample/runAnalyticalSolution.m @@ -26,7 +26,7 @@ %% Solution ifdiff integrator = @ode45; odeoptionsrhs_test = odeset( 'AbsTol', 1e-14,'RelTol', 1e-12); -datahandle = prepareDatahandleForIntegration('canonicalExampleRHS', 'integrator', func2str(integrator), 'options', odeoptionsrhs_test); +datahandle = prepareDatahandleForIntegration('rhsCanonicalExample', 'integrator', integrator, 'options', odeoptionsrhs_test); sol = solveODE(datahandle, tspan, y0, p); %% Plot analytical solution vs. ifdiff @@ -44,7 +44,7 @@ %set(gca,'XTick',0:5:20); %% Plot analytical solution vs ode45 -canonicalExampleRHS_ode45 = @(t,y) canonicalExampleRHS(t,y,p); +canonicalExampleRHS_ode45 = @(t,y) rhsCanonicalExample(t,y,p); sol_ode45 = ode45(canonicalExampleRHS_ode45, tspan, y0); figure diff --git a/toolbox/examples/canonicalExample/canonicalExampleRHS_test.m b/toolbox/examples/canonicalExample/runCanonicalExample.m similarity index 92% rename from toolbox/examples/canonicalExample/canonicalExampleRHS_test.m rename to toolbox/examples/canonicalExample/runCanonicalExample.m index 3771a592..fde1e5d6 100644 --- a/toolbox/examples/canonicalExample/canonicalExampleRHS_test.m +++ b/toolbox/examples/canonicalExample/runCanonicalExample.m @@ -4,13 +4,19 @@ timeinterval = [t0,tf]; initstates = [1 0 ]; p = 5.437; +odeoptions = odeset( 'AbsTol', 1e-20, 'RelTol', 1e-12, 'MaxStep', 1000); + +%% +fprintf('Integration with %s...\n ', func2str(integrator)); +th = tic(); +sol_matlab = integrator(@(t,x) rhsCanonicalExample(t,x,p), timeinterval, initstates, odeoptions); +time_matlab = toc(th); fprintf('Took %g seconds\n', time_matlab); %% fprintf('Preprocessing...\n '); -odeoptions = odeset( 'AbsTol', 1e-20, 'RelTol', 1e-12); -filename = 'canonicalExampleRHS'; +filename = 'rhsCanonicalExample'; th = tic(); -dhandle = prepareDatahandleForIntegration(filename, 'integrator', func2str(integrator), 'options', odeoptions); +dhandle = prepareDatahandleForIntegration(filename, 'integrator', integrator, 'options', odeoptions); time_prepare = toc(th); fprintf('Took %g seconds\n', time_prepare); %% @@ -19,13 +25,6 @@ sol_ifdiff = solveODE(dhandle, timeinterval, initstates, p); time_ifdiff = toc(th); fprintf('Took %g seconds\n', time_ifdiff); -%% -fprintf('Integration with %s...\n ', func2str(integrator)); -th = tic(); -sol_matlab = integrator(@(t,x) canonicalExampleRHS(t,x,p), timeinterval, initstates, odeoptions); -time_matlab = toc(th); fprintf('Took %g seconds\n', time_matlab); - - %% % do explicit euler integration fprintf('AccurateEuler integration...\n '); @@ -38,7 +37,7 @@ dt = (tf-t0) / N_euler; % time increment for k = 1:N_euler T(k+1) = t0 + k*dt; - X(:,k+1) = X(:,k) + dt * canonicalExampleRHS(T(k), X(:,k), p); + X(:,k+1) = X(:,k) + dt * rhsCanonicalExample(T(k), X(:,k), p); end time_euler = toc(th); fprintf('Took %g seconds\n', time_euler); skipper = floor(N_euler/1000); diff --git a/toolbox/examples/canonicalExample/sensitivities_CanonicalExample.m b/toolbox/examples/canonicalExample/sensCanonicalExample.m similarity index 95% rename from toolbox/examples/canonicalExample/sensitivities_CanonicalExample.m rename to toolbox/examples/canonicalExample/sensCanonicalExample.m index d8153004..ef874e85 100644 --- a/toolbox/examples/canonicalExample/sensitivities_CanonicalExample.m +++ b/toolbox/examples/canonicalExample/sensCanonicalExample.m @@ -2,7 +2,7 @@ integrator = @ode45; odeoptionsrhs_test = odeset( 'AbsTol', 1e-14,'RelTol', 1e-12); -datahandle = prepareDatahandleForIntegration('canonicalExampleRHS', 'integrator', func2str(integrator), 'options', odeoptionsrhs_test); +datahandle = prepareDatahandleForIntegration('rhsCanonicalExample', 'integrator', integrator, 'options', odeoptionsrhs_test); tspan = [0 20]; initialvalues = [1;0]; @@ -40,7 +40,7 @@ t_plot = 0:0.01:20; %sensitivities_END_plot = sensitivities_function_ENDpiecewise(t_plot); -sensitivities_END_plot = sensitivities_function_VDE(t_plot); +sensitivities_END_plot = sensitivities_function_VDE(t_plot); %correct? sensdata_y11 = arrayfun( @(x) x.Gy(1,1), sensitivities_END_plot); sensdata_y12 = arrayfun( @(x) x.Gy(1,2), sensitivities_END_plot); diff --git a/toolbox/examples/coulomb/coulomb_rhs.m b/toolbox/examples/coulomb/rhsCoulomb.m similarity index 96% rename from toolbox/examples/coulomb/coulomb_rhs.m rename to toolbox/examples/coulomb/rhsCoulomb.m index 8f1dda00..b36d2d6f 100644 --- a/toolbox/examples/coulomb/coulomb_rhs.m +++ b/toolbox/examples/coulomb/rhsCoulomb.m @@ -1,4 +1,4 @@ -function dx = coulomb_rhs(t,x,U) +function dx = rhsCoulomb(t,x,U) % INPUT: t - time % x - state % U - function handle of one variable (time) [voltage] diff --git a/toolbox/examples/coulomb/coulomb_test.m b/toolbox/examples/coulomb/runCoulomb.m similarity index 93% rename from toolbox/examples/coulomb/coulomb_test.m rename to toolbox/examples/coulomb/runCoulomb.m index 83814dc4..85cf5bf1 100644 --- a/toolbox/examples/coulomb/coulomb_test.m +++ b/toolbox/examples/coulomb/runCoulomb.m @@ -32,7 +32,7 @@ dt = (times.tf - times.t0) / Neuler; % time increment for k = 1:Neuler T(k+1) = times.t0 + k*dt; - X(:,k+1) = X(:,k) + dt * coulomb_rhs(T(k), X(:,k), Ufun); + X(:,k+1) = X(:,k) + dt * rhsCoulomb(T(k), X(:,k), Ufun); %if k > 10; break; end end fprintf('Euler integration took %g seconds.\n', toc()); @@ -43,7 +43,7 @@ % try with ode-solver odesolver = @ode15s; tic -sol = odesolver(@(t,x) coulomb_rhs(t,x,@U_umax3), [times.t0, times.tf], x0); +sol = odesolver(@(t,x) rhsCoulomb(t,x,@U_umax3), [times.t0, times.tf], x0); fprintf('%s integration took %g seconds.\n', func2str(odesolver), toc()); Tode = Teuler; Xode = deval(sol, Tode); diff --git a/toolbox/examples/daeExample/daeExample_main.m b/toolbox/examples/daeExample/daeExample_main.m deleted file mode 100644 index c840fa1f..00000000 --- a/toolbox/examples/daeExample/daeExample_main.m +++ /dev/null @@ -1,26 +0,0 @@ -% DAE Example Main -% for more details, see daeExample_README.md - -%% Setup and integration -integrator = @ode15s; -x0 = [1; -1]; -tspan = [0 5]; -M = [1 0; 0 0]; -p = -0.2; - -opts_ifdiff = odeset('Mass', M, 'MassSingular', 'yes', 'AbsTol', 1e-9,'RelTol', 1e-6); -opts_plain = odeset('Mass', M, 'MassSingular', 'yes', 'AbsTol', 1e-9, 'RelTol', 1e-6); - -datahandle = prepareDatahandleForIntegration('daeExampleRHS', 'integrator', integrator, 'options', opts_ifdiff); -sol_ifdiff = solveODE(datahandle, tspan, x0, p); -sol_plain = integrator(@(t, x) daeExampleRHS(t, x, p), tspan, x0, opts_plain); - -%% Plots -clf; -fig1 = figure(01); -hold on; -IFDIFF_plot_1 = plot(sol_ifdiff.x, sol_ifdiff.y, 'ro--', 'DisplayName', 'IFDIFF'); -Plain_plot_1 = plot(sol_plain.x, sol_plain.y, 'ko-', 'DisplayName', 'plain ode15s'); -Switch_plot = xline(sol_ifdiff.switches, 'b', 'LineWidth', 1.0, 'DisplayName', 'Switch'); -legend([Plain_plot_1(1), IFDIFF_plot_1(1), Switch_plot]); -hold off; diff --git a/toolbox/examples/daeExample/daeExample_README.md b/toolbox/examples/daeExamples/acadamicDaeExample/README_DaeExample.md similarity index 98% rename from toolbox/examples/daeExample/daeExample_README.md rename to toolbox/examples/daeExamples/acadamicDaeExample/README_DaeExample.md index 07c980ef..4afb48b7 100644 --- a/toolbox/examples/daeExample/daeExample_README.md +++ b/toolbox/examples/daeExamples/acadamicDaeExample/README_DaeExample.md @@ -123,12 +123,12 @@ legend([Plain_plot_1(1), IFDIFF_plot_1(1), Switch_plot]); hold off ``` -![](plots_daeExample/plot1.png) +![](plotsDaeExample/plot1.png) We notice that the integrator strategy results in small steps here for the plain solver as well as IFDIFF. This is standard behavior for `ode15s` which is a multi-step method; it is not a defect caused by improper treatment of switching events. However, if we take a closer look, we see that the integration with IFDIFF is accurate around the switching point. -![](plots_daeExample/plot1_close.png) +![](plotsDaeExample/plot1_close.png) ## Additional Content diff --git a/toolbox/examples/daeExample/plots_daeExample/plot1.png b/toolbox/examples/daeExamples/acadamicDaeExample/plotsDaeExample/plot1.png similarity index 100% rename from toolbox/examples/daeExample/plots_daeExample/plot1.png rename to toolbox/examples/daeExamples/acadamicDaeExample/plotsDaeExample/plot1.png diff --git a/toolbox/examples/daeExample/plots_daeExample/plot1_close.png b/toolbox/examples/daeExamples/acadamicDaeExample/plotsDaeExample/plot1_close.png similarity index 100% rename from toolbox/examples/daeExample/plots_daeExample/plot1_close.png rename to toolbox/examples/daeExamples/acadamicDaeExample/plotsDaeExample/plot1_close.png diff --git a/toolbox/examples/daeExample/daeExampleRHS.m b/toolbox/examples/daeExamples/acadamicDaeExample/rhsDaeExample.m similarity index 54% rename from toolbox/examples/daeExample/daeExampleRHS.m rename to toolbox/examples/daeExamples/acadamicDaeExample/rhsDaeExample.m index a13027ed..4b294305 100644 --- a/toolbox/examples/daeExample/daeExampleRHS.m +++ b/toolbox/examples/daeExamples/acadamicDaeExample/rhsDaeExample.m @@ -1,13 +1,16 @@ -function f = daeExampleRHS(~, x, p) - +function f = rhsDaeExample(~, x, p) + % f = zeros(2,1); + % algebraic constraint - z = x(1) + x(2); + z = x(1) + x(2) + x(2)^3 - p; f(2) = z; - % differential variables + + % differential equation if x(2) < p f(1) = x(2); else f(1) = 0; end -end + +end \ No newline at end of file diff --git a/toolbox/examples/daeExamples/acadamicDaeExample/runDaeExample.m b/toolbox/examples/daeExamples/acadamicDaeExample/runDaeExample.m new file mode 100644 index 00000000..dda38ed0 --- /dev/null +++ b/toolbox/examples/daeExamples/acadamicDaeExample/runDaeExample.m @@ -0,0 +1,34 @@ +% DAE Example Main +% for more details, see [new README missing] + +%% Setup and integration +integrator = @ode15s; +p = -0.4; % -1 < p < 0 +x0 = [1; -1+p]; % 1+p to guarantee consistent initial cond. +t0 = 0; +tf = 5; +tspan = [t0 tf]; +M = [1 0; 0 0]; + +opts_ifdiff = odeset('Mass', M, 'MassSingular', 'yes', 'AbsTol', 1e-8,'RelTol', 1e-6); +opts_plain = odeset('Mass', M, 'MassSingular', 'yes', 'AbsTol', 1e-8, 'RelTol', 1e-6); + +sol_plain = integrator(@(t, x) rhsDaeExample(t, x, p), tspan, x0, opts_plain); + +datahandle = prepareDatahandleForIntegration('rhsDaeExample', 'integrator', integrator, 'options', opts_ifdiff); +sol_ifdiff = solveODE(datahandle, tspan, x0, p); + +%% Expl. Euler comparision +h = 1e-4; +sol_euler = explicitEulerDAE(@(t,x,p) rhsDaeExample(t,x,p), 1, tspan, x0, p, h); + +%% Plots +clf; +fig2 = figure(01); +hold on; +Euler_plot = plot(sol_euler.x, sol_euler.y, 'g*--', 'DisplayName', 'Euler'); +IFDIFF_plot = plot(sol_ifdiff.x, sol_ifdiff.y, 'ro--', 'DisplayName', 'IFDIFF'); +Plain_plot = plot(sol_plain.x, sol_plain.y, 'ko-', 'DisplayName', 'plain ode15s'); +Switch_plot = xline(sol_ifdiff.switches, 'b', 'LineWidth', 1.0, 'DisplayName', 'Switch'); +legend([Plain_plot(1), IFDIFF_plot(1), Euler_plot(1), Switch_plot ]); +hold off; diff --git a/toolbox/examples/daeExamples/daeExampleWorkshop/README_DAE_Workshop.md b/toolbox/examples/daeExamples/daeExampleWorkshop/README_DAE_Workshop.md new file mode 100644 index 00000000..ef3fe8f0 --- /dev/null +++ b/toolbox/examples/daeExamples/daeExampleWorkshop/README_DAE_Workshop.md @@ -0,0 +1,20 @@ +# README Differential Algebraic Equation Example (from Workshop 2026) + +The mathematical formulation for ´rhsDaeExampleWorkshop.m´ +Let $t \in [0,3]$, $\text{x} = (x_1, x_2)^T$: + +$$ + (D) \quad + \begin{cases} + \dot{x}_ 1 = f(t,x,p) = + \begin{cases} + x_2, & \text{if } x_2 < p \\ + 0, & \text{if } x_2 \geq p + \end{cases} \\ + g(t,x,p) = x_1 + x_2 + 10(x_1^5 + x_2^3) = 0 \\ + \text{x}_0 = (1.0, -1.0)^T, p = -0.3 + \end{cases} +$$ + +On a larger scale, ´ode15s´ solution looks correct, but it fails around the switch at $t_s \approx 1.068$. When we zoom in, we see that ´ode15s´ does not notice the switch but instead fits a polynomial where the dynamics should change to constant immediately. +Since ´ode15s´ assumes a continuously differentiable solution trajectory and does not feature switching point detection, it is the wrong means of solving this problem. diff --git a/toolbox/examples/daeExamples/daeExampleWorkshop/rhsDaeExampleWorkshop.m b/toolbox/examples/daeExamples/daeExampleWorkshop/rhsDaeExampleWorkshop.m new file mode 100644 index 00000000..58055e61 --- /dev/null +++ b/toolbox/examples/daeExamples/daeExampleWorkshop/rhsDaeExampleWorkshop.m @@ -0,0 +1,16 @@ +function f = rhsDaeExampleWorkshop(~, x, p) + + f = zeros(2,1); + + % algebraic constraint + z = x(1) + x(2) + 10*(x(1)^5 + x(2)^3); + f(2) = z; + + % differential equation + if x(2) < p + f(1) = x(2); + else + f(1) = 0; + end + +end diff --git a/toolbox/examples/daeExamples/daeExampleWorkshop/runDaeExampleWorkshop.m b/toolbox/examples/daeExamples/daeExampleWorkshop/runDaeExampleWorkshop.m new file mode 100644 index 00000000..c59ae7cd --- /dev/null +++ b/toolbox/examples/daeExamples/daeExampleWorkshop/runDaeExampleWorkshop.m @@ -0,0 +1,60 @@ +% Differential Algebraic Equation Example +% from Workshop 2026 (more info, see README_DAE_Workshop.md) + +%% Setup +integrator = @ode15s; +x0 = [1; -1]; +t0 = 0; +tf = 3; +tspan = [t0 tf]; +p = -0.3; + +% mass matrix +M = [1 0; 0 0]; +opts = odeset('Mass', M, 'MassSingular', 'yes', 'AbsTol', 1e-9,'RelTol', 1e-5); + + +%% Solution with IFDIFF +filename = 'rhsDaeExampleWorkshop'; +datahandle = prepareDatahandleForIntegration(filename, 'integrator', integrator, 'options', opts); +sol_ifdiff = solveODE(datahandle, tspan, x0, p); + +%% Solution with ode15s +sol_plain = integrator(@(t, x) rhsDaeExampleWorkshop(t, x, p), tspan, x0, opts); + + +%% Plots + +t_eval = t0:1e-6:tf; +y_ifdiff = deval(sol_ifdiff, t_eval); +y_plain = deval(sol_plain, t_eval); + +% plot colors +c_ifdiff = [0.85, 0.33, 0.1]; +c_plain = [0.0, 0.45, 0.74]; + + +% large scale plot +% on a larger scale both solutions look correct +fig1 = figure(1); +clf(fig1, "reset"); +hold("on"); +IFDIFF_plot = plot(t_eval, y_ifdiff, '-', 'Color', c_ifdiff, 'LineWidth', 1.2, 'DisplayName', 'IFDIFF'); +Plain_plot = plot(t_eval, y_plain, '--', 'Color', c_plain, 'DisplayName', 'plain ode15s'); +legend([Plain_plot(1), IFDIFF_plot(1)]); +grid("on"); +hold("off"); + +% zoom in at switching point +% ode15s interpolates around the switch whereas ifdiff detects and handles +% it correctly +fig2 = figure(2); +clf(fig2, "reset"); +hold("on"); +IFDIFF_plot = plot(t_eval, y_ifdiff, '-', 'Color', c_ifdiff, 'DisplayName', 'IFDIFF'); +Plain_plot = plot(t_eval, y_plain, '-', 'Color', c_plain, 'DisplayName', 'plain ode15s'); +xlim([1.0664 1.0702]) +ylim([0.4271 0.4277]) +legend([Plain_plot(1), IFDIFF_plot(1)]); +grid("on"); +hold("off"); diff --git a/toolbox/examples/daeExamples/explicitEulerDAE.m b/toolbox/examples/daeExamples/explicitEulerDAE.m new file mode 100644 index 00000000..61e0ad18 --- /dev/null +++ b/toolbox/examples/daeExamples/explicitEulerDAE.m @@ -0,0 +1,75 @@ +function sol = explicitEulerDAE(rhs, DifferentialVars, tspan, x0, p, h, newtonOpts) +% explicitEulerDAE - Explicit Euler solver for semi-explicit, DAEs of +% differential index 1 +% +% Solves: M * x' = f(t,x) by explicit Euler method +% INPUT: +% rhs - function handle: f(t, x, p) +% tspan - [t0 tf] +% x0 - initial state +% p - parameters +% h - step size +% newtonOpts - struct with field +% DifferentialVars - number of differential variables +% OPTIONS: +% maxIter (default 10) +% tol (default 1e-10) +% OUTPUT: +% sol.x - time vector +% sol.y - solution matrix + + % default options + if nargin < 7 + newtonOpts.maxIter = 10; + newtonOpts.tol = 1e-10; + end + + t = tspan(1):h:tspan(2); + N = length(t); + + nx = length(x0); + x = zeros(nx, N); + x(:,1) = x0; + + % calculate dimensions + nd = DifferentialVars; % number of differential variables + na = nx - nd; % number of algebraic variables + + for n = 1:N-1 + f_val = rhs(t(n), x(:,n), p); + x(1:nd, n+1) = x(1:nd, n) + h * f_val(1:nd); + z_guess = x(nd+1:end, n); + + for k = 1:newtonOpts.maxIter + x_trial = [x(1:nd, n+1); z_guess]; + F = rhs(t(n+1), x_trial, p); + g = F(nd+1:end); + J = zeros(na, na); + eps_fd = 1e-8; + + for i = 1:na + dz = zeros(na,1); + dz(i) = eps_fd; + + x_eps = [x(1:nd, n+1); z_guess + dz]; + F_eps = rhs(t(n+1), x_eps, p); + J(:,i) = (F_eps(nd+1:end) - g) / eps_fd; + end + + % Newton step (differential variables) + update = -J \ g; + z_new = z_guess + update; + + if norm(update) < newtonOpts.tol + break; + end + + z_guess = z_new; + end + + x(nd+1:end, n+1) = z_guess; + end + + sol.x = t; + sol.y = x; +end \ No newline at end of file diff --git a/toolbox/examples/rlcExample/plots_rlc/rlc_plot.png b/toolbox/examples/daeExamples/rlcExample/plotsRLC/rlc_plot.png similarity index 100% rename from toolbox/examples/rlcExample/plots_rlc/rlc_plot.png rename to toolbox/examples/daeExamples/rlcExample/plotsRLC/rlc_plot.png diff --git a/toolbox/examples/rlcExample/plots_rlc/rlc_plot_close.png b/toolbox/examples/daeExamples/rlcExample/plotsRLC/rlc_plot_close.png similarity index 100% rename from toolbox/examples/rlcExample/plots_rlc/rlc_plot_close.png rename to toolbox/examples/daeExamples/rlcExample/plotsRLC/rlc_plot_close.png diff --git a/toolbox/examples/rlcExample/rlcRHS.m b/toolbox/examples/daeExamples/rlcExample/rhsRLC.m similarity index 91% rename from toolbox/examples/rlcExample/rlcRHS.m rename to toolbox/examples/daeExamples/rlcExample/rhsRLC.m index d601beeb..e76977b5 100644 --- a/toolbox/examples/rlcExample/rlcRHS.m +++ b/toolbox/examples/daeExamples/rlcExample/rhsRLC.m @@ -1,4 +1,4 @@ -function dx = rlcRHS(~,x,p) +function dx = rhsRLC(~,x,p) dx = zeros(3,1); L = p(1); R1 = p(2); R2 = p(3); C = p(4); Vs = p(5); Vth = p(6); diff --git a/toolbox/examples/rlcExample/rlcExample_README.md b/toolbox/examples/daeExamples/rlcExample/rlcExample_README.md similarity index 100% rename from toolbox/examples/rlcExample/rlcExample_README.md rename to toolbox/examples/daeExamples/rlcExample/rlcExample_README.md diff --git a/toolbox/examples/rlcExample/rlc_main.m b/toolbox/examples/daeExamples/rlcExample/runRLC.m similarity index 58% rename from toolbox/examples/rlcExample/rlc_main.m rename to toolbox/examples/daeExamples/rlcExample/runRLC.m index 0e5e5f20..703ab6a5 100644 --- a/toolbox/examples/rlcExample/rlc_main.m +++ b/toolbox/examples/daeExamples/rlcExample/runRLC.m @@ -18,10 +18,15 @@ opts_ifdiff = odeset('Mass', M,'AbsTol', 1e-8, 'RelTol', 1e-5); opts_plain = odeset('Mass', M,'AbsTol', 1e-8, 'RelTol', 1e-5); -datahandle = prepareDatahandleForIntegration('rlcRHS', 'integrator', integrator, 'options', opts_ifdiff); +datahandle = prepareDatahandleForIntegration('rhsRLC', 'integrator', integrator, 'options', opts_ifdiff); sol_ifdiff = solveODE(datahandle, tspan, x0, p); -sol_plain = integrator(@(t, x) rlcRHS(t, x, p), tspan, x0, opts_plain); +sol_plain = integrator(@(t, x) rhsRLC(t, x, p), tspan, x0, opts_plain); + + +%% Explicit Euler (for comparison) +h = 1e-4; +sol_euler = explicitEulerDAE(@(t,x,p) rhsRLC(t,x,p), 2, tspan, x0, p, h); %% Plots clf; @@ -29,18 +34,20 @@ subplot(2,1,1); hold on; +Euler_plot_1 = plot(sol_euler.x, sol_euler.y(1,:), 'go--', 'DisplayName', 'Explicit Euler'); Plot_ifdiff_1 = plot(sol_ifdiff.x, sol_ifdiff.y(1,:), 'ro--', 'DisplayName', 'IFDIFF'); Plot_plain_1 = plot(sol_plain.x, sol_plain.y(1,:), 'k.-', 'DisplayName', 'plain ode15s'); -%Switch_plot_1 = xline(sol_ifdiff.switches, 'b', 'LineWidth', 1.0, 'DisplayName', 'Switch'); +Switch_plot_1 = xline(sol_ifdiff.switches, 'b', 'LineWidth', 1.0, 'DisplayName', 'Switch'); ylabel('i_L (A)'); xlabel('Time (s)'); legend(); hold off; subplot(2,1,2); -hold on -Plot_ifdiff_2 = plot(sol_ifdiff.x, sol_ifdiff.y(2,:), 'ro--', 'DisplayName', 'IFDIFF' ); -Plot_plain_2 = plot(sol_plain.x, sol_plain.y(2,:), 'k.-', 'DisplayName', 'plain ode15s'); +hold on; +Euler_plot_2 = plot(sol_euler.x, sol_euler.y(2,:), 'go--', 'DisplayName', 'Explicit Euler'); +Plot_ifdiff_2 = plot(sol_ifdiff.x, sol_ifdiff.y(2,:), 'ro--', 'DisplayName', 'IFDIFF' ); +Plot_plain_2 = plot(sol_plain.x, sol_plain.y(2,:), 'k.-', 'DisplayName', 'plain ode15s'); ylabel('i_C (A) = i_L (A) (via constraint)'); Switch_plot_2 = xline(sol_ifdiff.switches, 'b', 'LineWidth', 1.0, 'DisplayName', 'Switch'); xlabel('Time (s)'); diff --git a/toolbox/examples/funfun/funinif_ext.m b/toolbox/examples/funfun/funinif_ext.m deleted file mode 100644 index 51f33346..00000000 --- a/toolbox/examples/funfun/funinif_ext.m +++ /dev/null @@ -1,3 +0,0 @@ -function val = funinif_ext(t,x) - val = exp(1)*(t+0)*x; -end \ No newline at end of file diff --git a/toolbox/examples/funfun/funinif_rhs.m b/toolbox/examples/funfun/funinif_rhs.m deleted file mode 100644 index 93203a9c..00000000 --- a/toolbox/examples/funfun/funinif_rhs.m +++ /dev/null @@ -1,7 +0,0 @@ -function val = funinif_rhs(t,y,~) - if t > funinif_ext(t,y) - val = 1; - else - val = 1.1; - end -end \ No newline at end of file diff --git a/toolbox/examples/funfun/funinif_test.m b/toolbox/examples/funfun/funinif_test.m deleted file mode 100644 index 3d4e3693..00000000 --- a/toolbox/examples/funfun/funinif_test.m +++ /dev/null @@ -1,40 +0,0 @@ -% Test if function call in if condition is allowed - -% initialize -x0 = 1; -t0 = 0; -tf = 9; -tspan = [t0 tf]; -T = t0:((tf-t0)/100):tf; - -% setup -odesolver = @ode45; -options = odeset('AbsTol', 1e-15, 'RelTol', 1e-3); -% options = odeset('AbsTol', 1e-20, 'RelTol', 1e-8, 'MaxStep', 1e-3); -rhs_name = 'funinif_rhs'; -% rhs_name = 'funinif_rhs'; - -% try with ode-solver -tic -solode = odesolver(@(t,x) funinif_rhs(t,x), tspan, x0, options); -fprintf('%s integration took %g seconds.\n', func2str(odesolver), toc()); -X = deval(solode, T); - -% ifdiff -tic -hdl = prepareDatahandleForIntegration(rhs_name, 'integrator', func2str(odesolver), 'options', options); -fprintf('Ifdiff handle preparation took %g seconds.\n', toc()); -tic -[solif, dataif] = solveODE(hdl, tspan, x0, []); -fprintf('Ifdiff integration took %g seconds.\n', toc()); -XX = deval(solif, T); - -% plot solutions -figure(313); clf; -subplot(2,1,1); plot(T, XX, 'go-'); hold on; -subplot(2,1,1); plot(T, X , 'b.-'); -legend('ifdiff',func2str(odesolver),'Location','Best'); -subplot(2,1,2); plot(T, X-XX, 'k.-'); -legend('difference','Location','Best'); - - diff --git a/toolbox/examples/liveExamples/CanonicalSensitivities.mlx b/toolbox/examples/liveExamples/CanonicalSensitivities.mlx deleted file mode 100644 index d4486ded..00000000 Binary files a/toolbox/examples/liveExamples/CanonicalSensitivities.mlx and /dev/null differ diff --git a/toolbox/examples/liveExamples/FilippovLive.mlx b/toolbox/examples/liveExamples/FilippovLive.mlx index 1b2878b2..92794a2a 100644 Binary files a/toolbox/examples/liveExamples/FilippovLive.mlx and b/toolbox/examples/liveExamples/FilippovLive.mlx differ diff --git a/toolbox/examples/liveExamples/StateJumpExample.mlx b/toolbox/examples/liveExamples/StateJumpExample.mlx index 5685e8a3..75ccb5e1 100644 Binary files a/toolbox/examples/liveExamples/StateJumpExample.mlx and b/toolbox/examples/liveExamples/StateJumpExample.mlx differ diff --git a/toolbox/examples/liveExamples/WhiteCabbageLive.mlx b/toolbox/examples/liveExamples/WhiteCabbageLive.mlx index 7b5ba968..3e409cb2 100644 Binary files a/toolbox/examples/liveExamples/WhiteCabbageLive.mlx and b/toolbox/examples/liveExamples/WhiteCabbageLive.mlx differ diff --git a/toolbox/examples/liveExamples/canonicalExRHS_test_live.mlx b/toolbox/examples/liveExamples/canonicalExRHS_test_live.mlx deleted file mode 100644 index 23d2bfb3..00000000 Binary files a/toolbox/examples/liveExamples/canonicalExRHS_test_live.mlx and /dev/null differ diff --git a/toolbox/examples/liveExamples/helperFunctions/plot_canon_sens_p.m b/toolbox/examples/liveExamples/helperFunctions/plotCanonSensParameter.m similarity index 94% rename from toolbox/examples/liveExamples/helperFunctions/plot_canon_sens_p.m rename to toolbox/examples/liveExamples/helperFunctions/plotCanonSensParameter.m index ecf9fd19..a83bcff5 100644 --- a/toolbox/examples/liveExamples/helperFunctions/plot_canon_sens_p.m +++ b/toolbox/examples/liveExamples/helperFunctions/plotCanonSensParameter.m @@ -1,4 +1,4 @@ -function plot_canon_sens_p(t_plot, sensitivities_END_plot, sensitivities_END_plot2) +function plotCanonSensParameter(t_plot, sensitivities_END_plot, sensitivities_END_plot2) % just for very specific example (only used in one live script) % arguments are as follows: %t_plot = 0:0.01:20; diff --git a/toolbox/examples/liveExamples/helperFunctions/plot_canonical_comparison.m b/toolbox/examples/liveExamples/helperFunctions/plotCanonicalComparison.m similarity index 86% rename from toolbox/examples/liveExamples/helperFunctions/plot_canonical_comparison.m rename to toolbox/examples/liveExamples/helperFunctions/plotCanonicalComparison.m index b6767c4e..7bb25701 100644 --- a/toolbox/examples/liveExamples/helperFunctions/plot_canonical_comparison.m +++ b/toolbox/examples/liveExamples/helperFunctions/plotCanonicalComparison.m @@ -1,4 +1,4 @@ -function plot_canonical_comparison(fignum, tf, T, Y_ifdiff, Y_euler, Y_matlabsolver, sol_matlab, sol_ifdiff, ymax) +function plotCanonicalComparison(fignum, tf, T, Y_ifdiff, Y_euler, Y_matlabsolver, sol_matlab, sol_ifdiff, ymax) figure(fignum); clf('reset'); hold on; diff --git a/toolbox/examples/liveExamples/helperFunctions/plot_comparison_subplots.m b/toolbox/examples/liveExamples/helperFunctions/plotComparisonSubplots.m similarity index 90% rename from toolbox/examples/liveExamples/helperFunctions/plot_comparison_subplots.m rename to toolbox/examples/liveExamples/helperFunctions/plotComparisonSubplots.m index 34a17fab..5497f7d2 100644 --- a/toolbox/examples/liveExamples/helperFunctions/plot_comparison_subplots.m +++ b/toolbox/examples/liveExamples/helperFunctions/plotComparisonSubplots.m @@ -1,4 +1,4 @@ -function plot_comparison_subplots(fignum, tf, t0, T, Y_ifdiff, Y_euler, Y_matlabsolver, sol_ifdiff, integrator) +function plotComparisonSubplots(fignum, tf, t0, T, Y_ifdiff, Y_euler, Y_matlabsolver, sol_ifdiff, integrator) Ydiff_Euler = Y_euler - Y_ifdiff; Ydiff_Matlab = Y_matlabsolver - Y_ifdiff; fignum = fignum + 1; diff --git a/toolbox/examples/liveExamples/helperFunctions/bounceballRHS_cp.m b/toolbox/examples/liveExamples/helperFunctions/rhsBounceballCopy.m similarity index 82% rename from toolbox/examples/liveExamples/helperFunctions/bounceballRHS_cp.m rename to toolbox/examples/liveExamples/helperFunctions/rhsBounceballCopy.m index 1d233a1f..ade9ef71 100644 --- a/toolbox/examples/liveExamples/helperFunctions/bounceballRHS_cp.m +++ b/toolbox/examples/liveExamples/helperFunctions/rhsBounceballCopy.m @@ -1,4 +1,4 @@ -function dx = bounceballRHS_cp(~, x, p) +function dx = rhsBounceballCopy(~, x, p) dx = [x(2); -p(1)]; if ifdiff_jumpif(x(1), -1) deltaH = -x(1) + eps(1)*(1/p(1)) * p(2)^2*x(2)^2; diff --git a/toolbox/examples/liveExamples/helperFunctions/run_accurate_euler.m b/toolbox/examples/liveExamples/helperFunctions/runAccurateEuler.m similarity index 77% rename from toolbox/examples/liveExamples/helperFunctions/run_accurate_euler.m rename to toolbox/examples/liveExamples/helperFunctions/runAccurateEuler.m index 7e8c04ee..14c1c890 100644 --- a/toolbox/examples/liveExamples/helperFunctions/run_accurate_euler.m +++ b/toolbox/examples/liveExamples/helperFunctions/runAccurateEuler.m @@ -1,4 +1,4 @@ -function [T, Y_euler, Y_matlabsolver, Y_ifdiff] = run_accurate_euler(N_euler, t0, tf, initstates, p, canonicalExampleRHS, sol_matlab, sol_ifdiff) +function [T, Y_euler, Y_matlabsolver, Y_ifdiff] = runAccurateEuler(N_euler, t0, tf, initstates, p, canonicalExampleRHS, sol_matlab, sol_ifdiff) X = zeros(2, N_euler+1); X(:,1) = initstates; T = zeros(1, N_euler+1); diff --git a/toolbox/examples/liveExamples/helperFunctions/simulateBouncingBall.m b/toolbox/examples/liveExamples/helperFunctions/simulateBouncingBall.m index 06d5b47d..e1ec5720 100644 --- a/toolbox/examples/liveExamples/helperFunctions/simulateBouncingBall.m +++ b/toolbox/examples/liveExamples/helperFunctions/simulateBouncingBall.m @@ -1,5 +1,5 @@ function simulateBouncingBall(t0, tEnd, x0, p, solver, options) - datahandle = prepareDatahandleForIntegration('bounceballRHS', 'solver', func2str(solver), 'options', options); + datahandle = prepareDatahandleForIntegration('rhsBounceball', 'solver', solver, 'options', options); sol = solveODE(datahandle, [t0 tEnd], x0, p); sensFun = generateSensitivityFunction(datahandle, sol, 'method', 'VDE', 'CalcGy', true, 'CalcGp', true); diff --git a/toolbox/examples/liveExamples/runCanonicalExampleLive.mlx b/toolbox/examples/liveExamples/runCanonicalExampleLive.mlx new file mode 100644 index 00000000..1f576895 Binary files /dev/null and b/toolbox/examples/liveExamples/runCanonicalExampleLive.mlx differ diff --git a/toolbox/examples/liveExamples/sensCanonicalExample.mlx b/toolbox/examples/liveExamples/sensCanonicalExample.mlx new file mode 100644 index 00000000..eadd59cd Binary files /dev/null and b/toolbox/examples/liveExamples/sensCanonicalExample.mlx differ diff --git a/toolbox/examples/parameterEstimationWhiteCabbage/HelpersWhiteCabbage/computeCI.m b/toolbox/examples/parameterEstimationWhiteCabbage/HelpersWhiteCabbage/computeCI.m new file mode 100644 index 00000000..cf733ca8 --- /dev/null +++ b/toolbox/examples/parameterEstimationWhiteCabbage/HelpersWhiteCabbage/computeCI.m @@ -0,0 +1,42 @@ +function CI = computeCI(p, res, J, alpha) +% CI = computeCI(p, res, J, varargin) +% computes the confidence confidence intervals in the parameter +% estimation problem in runCabbage.m +% INPUT: +% p - parameter vector (of optimized parameters) +% res - residual (from lsqnonlin solve) +% J - Jacobian +% varargin: alpha - parameter in (0,1) +% OUTPUT: +% CI - two-sided confidence intervals + +if nargin < 4 + alpha = 0.05; % defaultfor 95% CI +end + +np = numel(p); +nr = numel(res); +dof = nr - np; +res_variance = sum(res.^2)/dof; + +% approximation when a column is zero vector +index = find(max(abs(J)) == 0); +if ~isempty(index) + J(:,index) = 1e-8; +end + +% computation of parameter covariance matrix by QR factorization +[~, R] = qr(J, 0); +R_inv = R \ eye(size(R)); +Cov_p = res_variance * R_inv * transpose(R_inv); + +% critical t-value at 1-alpha/2 +t_crit = tInv(alpha, dof); + +% standard error +std_errors = sqrt(diag(Cov_p)); +margin = t_crit * std_errors; + +CI = [p - margin, p + margin]; + +end diff --git a/toolbox/examples/parameterEstimationWhiteCabbage/generateResidualFunction.m b/toolbox/examples/parameterEstimationWhiteCabbage/HelpersWhiteCabbage/generateResidualFunction.m similarity index 91% rename from toolbox/examples/parameterEstimationWhiteCabbage/generateResidualFunction.m rename to toolbox/examples/parameterEstimationWhiteCabbage/HelpersWhiteCabbage/generateResidualFunction.m index da12e424..e7775570 100644 --- a/toolbox/examples/parameterEstimationWhiteCabbage/generateResidualFunction.m +++ b/toolbox/examples/parameterEstimationWhiteCabbage/HelpersWhiteCabbage/generateResidualFunction.m @@ -7,7 +7,7 @@ integrator_solveODE = @ode45; odeoptionsrhs_test = odeset( 'AbsTol', 1e-14,'RelTol', 1e-6); - datahandle = prepareDatahandleForIntegration('whiteCabbageRHS', 'integrator', func2str(integrator_solveODE), 'options', odeoptionsrhs_test); + datahandle = prepareDatahandleForIntegration('rhsCabbage', 'integrator', integrator_solveODE, 'options', odeoptionsrhs_test); return diff --git a/toolbox/examples/parameterEstimationWhiteCabbage/getParams_Cabbage.m b/toolbox/examples/parameterEstimationWhiteCabbage/HelpersWhiteCabbage/getParamsCabbage.m similarity index 90% rename from toolbox/examples/parameterEstimationWhiteCabbage/getParams_Cabbage.m rename to toolbox/examples/parameterEstimationWhiteCabbage/HelpersWhiteCabbage/getParamsCabbage.m index 82f94cf8..76e3d821 100644 --- a/toolbox/examples/parameterEstimationWhiteCabbage/getParams_Cabbage.m +++ b/toolbox/examples/parameterEstimationWhiteCabbage/HelpersWhiteCabbage/getParamsCabbage.m @@ -1,4 +1,4 @@ -function parameters = getParams_Cabbage() +function parameters = getParamsCabbage() %Parameters for growth of white cabbage model parameters(1) = 34469.4; %a diff --git a/toolbox/examples/parameterEstimationWhiteCabbage/HelpersWhiteCabbage/tInv.m b/toolbox/examples/parameterEstimationWhiteCabbage/HelpersWhiteCabbage/tInv.m new file mode 100644 index 00000000..c992251b --- /dev/null +++ b/toolbox/examples/parameterEstimationWhiteCabbage/HelpersWhiteCabbage/tInv.m @@ -0,0 +1,47 @@ +function t_crit = tInv(alpha, dof) +% tcrit = tInv(alpha, dof) +% Computes the two-sided critical t-value for significance level alpha. +% INPUT: +% alpha - significance level in (0,1) (e.g. alpha = 0.05 for 95% CI) +% dof - degrees of freedom +% OUTPUT: +% t_crit - critical t-value + +% two sided confidence level (e.g. 0.975 if alpha = 0.05) +conf_level = 1 - 0.5*alpha; + +% input validation +if dof <= 0 + error('Degrees of freedom must be positive'); +end +if conf_level < 0 || conf_level > 1 + error('Confidence level must be in (0,1)'); +end + +% at mid-point +tol = 1e-12; +if abs(conf_level - 0.5) < tol + t_crit = 0; + return; +end + +% computation for negative t (due to point symmetry of CDF) +if conf_level < 0.5 + sign_t = -1; + conf_level_pos = 1 - conf_level; +else + sign_t = 1; + conf_level_pos = conf_level; +end + +% determine t_crit +% CDF(t) = 1 - 0.5*betainc(dof/(dof+t^2), dof/2, 0.5) +% Solve: CDF(t) = conf_level_pos + +objFunc = @(t) 1 - 0.5*betainc(dof/(dof + t^2), dof/2, 0.5) - conf_level_pos; + +options = optimoptions('fsolve', 'TolFun', 1e-10, 'Display', 'off'); +t_pos = fsolve(objFunc, 1, options); +t_crit = sign_t * t_pos; + +end \ No newline at end of file diff --git a/toolbox/examples/parameterEstimationWhiteCabbage/plot_sensitivities_cabbage.m b/toolbox/examples/parameterEstimationWhiteCabbage/plotSensitivitiesCabbage.m similarity index 68% rename from toolbox/examples/parameterEstimationWhiteCabbage/plot_sensitivities_cabbage.m rename to toolbox/examples/parameterEstimationWhiteCabbage/plotSensitivitiesCabbage.m index 62569868..55fee6c7 100644 --- a/toolbox/examples/parameterEstimationWhiteCabbage/plot_sensitivities_cabbage.m +++ b/toolbox/examples/parameterEstimationWhiteCabbage/plotSensitivitiesCabbage.m @@ -1,12 +1,22 @@ %% Solution ifdiff integrator = @ode45; +<<<<<<< HEAD +odeoptionsrhs_test = odeset( 'AbsTol', 1e-14,'RelTol', 1e-6); +datahandle = prepareDatahandleForIntegration('rhsCabbage', 'integrator', integrator, 'options', odeoptionsrhs_test); + + +tspan = [0 118]; +initialvalues = [2.48252;0;0]; +parameters_ODE = getParamsCabbage(); +======= odeoptionsrhs_test = odeset( 'AbsTol', 1e-14,'RelTol', 1e-12); datahandle = prepareDatahandleForIntegration('whiteCabbageRHS', 'integrator', func2str(integrator), 'options', odeoptionsrhs_test); tspan = [0 118]; initialvalues = [2.48252;0;0]; parameters_ODE = getParams_Cabbage(); +>>>>>>> 60ae785 (rename in White Cabbage and Predator Prey Filippov) sol = solveODE(datahandle, tspan, initialvalues, parameters_ODE); %% Precalculations for finite differences for sensitivities @@ -24,6 +34,45 @@ %% figure(1) +<<<<<<< HEAD +clf + +% Extract Gy matrices +GyCell = arrayfun(@(s) s.Gy, sensitivites, 'UniformOutput', false); +Gy = cat(3, GyCell{:}); + +labels = 'LSH'; +set(gcf, 'Position', [100 100 1100 850]); + +for row = 1:3 + for col = 1:3 + subplot(3,3,(row-1)*3 + col) + + plot(t_sens, squeeze(Gy(row,col,:)), '.', 'Color', [0.1 0.35 0.8], 'MarkerSize', 7); + + if col == 1 + ylim([0 500]) + else + ylim([-1 1.5]) + end + xlim([0 118]) + + xlabel('$t$', 'Interpreter', 'latex', 'FontSize', 13) + ylabel(sprintf('$\\partial y_{%s}(t)/\\partial y_{0,%s}$', labels(row), labels(col)), 'Interpreter', 'latex', 'FontSize', 13) + title(sprintf('$G_{y,%d%d}(t;t_0)$', row, col), 'Interpreter', 'latex', 'FontSize', 14, 'FontWeight', 'normal') + set(gca, ... + 'FontSize', 11, ... + 'LineWidth', 1, ... + 'Box', 'off', ... + 'TickDir', 'out', ... + 'TickLength', [0.015 0.015]) + + grid on + set(gca, 'GridAlpha', 0.15) + + end +end +======= subplot(3,3,1) hold on for i = 1:length(t_sens) @@ -140,10 +189,60 @@ title('G_{y,33}(t; t_0)') set(gca, 'FontSize', 22); set(gca, 'Box', 'off'); +>>>>>>> 60ae785 (rename in White Cabbage and Predator Prey Filippov) %% Plot sensitivities parameters sensitivities_function_p = generateSensitivityFunction(datahandle, sol, 'FDstep', FDstep, 'calcGy', false, 'method', method, 'p_typ', parameters_ODE); t_sens = 0:0.1:118; +<<<<<<< HEAD +ticid = tic; +sensitivites = sensitivities_function_p(t_sens); +toc(ticid); +%% +figure(2) +clf + +% Extract all Gp matrices +GpCell = arrayfun(@(s) s.Gp, sensitivites, 'UniformOutput', false); +Gp = cat(3, GpCell{:}); + +labels = 'LSH'; + +cols = [1 5 6]; + +ylims = [ + 0 0.02 + -20 150 + 0 4000 +]; + +paramLabels = {'a', 'r_S', 'r_H'}; +set(gcf, 'Position', [100 100 1100 850]); + +for row = 1:3 + for j = 1:3 + col = cols(j); + subplot(3,3,(row-1)*3 + j) + plot(t_sens, squeeze(Gp(row,col,:)), '.', 'Color', [0.1 0.35 0.8], 'MarkerSize', 7); + + xlim([0 118]) + ylim(ylims(j,:)) + + xlabel('$t$', 'Interpreter', 'latex', 'FontSize', 13) + ylabel(sprintf('$\\partial y_{%s}(t)/\\partial %s$', labels(row), paramLabels{j}), 'Interpreter', 'latex', 'FontSize', 13) + title(sprintf('$G_{p,%d%d}(t;t_0)$', row, col), 'Interpreter', 'latex', 'FontSize', 14, 'FontWeight', 'normal') + set(gca, ... + 'FontSize', 11, ... + 'LineWidth', 1, ... + 'Box', 'off', ... + 'TickDir', 'out', ... + 'TickLength', [0.015 0.015]) + + grid on + set(gca, 'GridAlpha', 0.15) + end +end +======= sensitivites = sensitivities_function_p(t_sens); %% @@ -264,3 +363,4 @@ title('G_{p,36}(t; t_0)') set(gca, 'FontSize', 22); set(gca, 'Box', 'off'); +>>>>>>> 60ae785 (rename in White Cabbage and Predator Prey Filippov) diff --git a/toolbox/examples/parameterEstimationWhiteCabbage/whiteCabbageRHS.m b/toolbox/examples/parameterEstimationWhiteCabbage/rhsCabbage.m similarity index 88% rename from toolbox/examples/parameterEstimationWhiteCabbage/whiteCabbageRHS.m rename to toolbox/examples/parameterEstimationWhiteCabbage/rhsCabbage.m index f90c5c38..42405e1c 100644 --- a/toolbox/examples/parameterEstimationWhiteCabbage/whiteCabbageRHS.m +++ b/toolbox/examples/parameterEstimationWhiteCabbage/rhsCabbage.m @@ -1,4 +1,4 @@ -function dy = whiteCabbageRHS(t, y, p) +function dy = rhsCabbage(t, y, p) %RHS for growth of white cabbage model dy = zeros(3,1); diff --git a/toolbox/examples/parameterEstimationWhiteCabbage/test_Cabbage.m b/toolbox/examples/parameterEstimationWhiteCabbage/runCabbage.m similarity index 81% rename from toolbox/examples/parameterEstimationWhiteCabbage/test_Cabbage.m rename to toolbox/examples/parameterEstimationWhiteCabbage/runCabbage.m index abaa956c..1bfc6d75 100644 --- a/toolbox/examples/parameterEstimationWhiteCabbage/test_Cabbage.m +++ b/toolbox/examples/parameterEstimationWhiteCabbage/runCabbage.m @@ -1,14 +1,13 @@ %% Solution ifdiff integrator = @ode45; -odeoptionsrhs_test = odeset( 'AbsTol', 1e-14,'RelTol', 1e-12); -datahandle = prepareDatahandleForIntegration('whiteCabbageRHS', 'integrator', func2str(integrator), 'options', odeoptionsrhs_test); +odeoptionsrhs = odeset( 'AbsTol', 1e-14,'RelTol', 1e-12); +datahandle = prepareDatahandleForIntegration('rhsCabbage', 'integrator', integrator, 'options', odeoptionsrhs); tspan = [0 118]; initialvalues = [2.48252;0;0]; -parameters_ODE = getParams_Cabbage(); +parameters_ODE = getParamsCabbage(); sol = solveODE(datahandle, tspan, initialvalues, parameters_ODE); -%s = rng; %% Precalculations for finite differences for sensitivities dim_y = size(sol.y, 1); @@ -19,7 +18,8 @@ t = 0:118; disturb = true; sigma = 5; -%rng(s); +seed = 72; +rng(seed); rand = randn(dim_y*length(t),1); if disturb measurements = reshape(deval(sol, t), [], 1) + sigma.*rand; @@ -29,9 +29,8 @@ %% Generation of residual function integrator_residual = @solveODE; -RHS = @whiteCabbageRHS; method = 'VDE'; -residual_function = generateResidualFunction(t, datahandle, sol, measurements, tspan, parameters_ODE, RHS, FDstep, integrator_residual, method); +residual_function = generateResidualFunction(t, datahandle, sol, measurements, tspan, parameters_ODE, @rhsCabbage, FDstep, integrator_residual, method); %% Parameter estimation options = optimoptions('lsqnonlin','SpecifyObjectiveGradient',true, 'Algorithm','levenberg-marquardt', 'Display', 'iter', 'DerivativeCheck', 'off', 'typicalX', [parameters_ODE;1;1;1]);%, 'TolX', 1e-14, 'TolFun', 1e-14); @@ -39,12 +38,14 @@ tic; [param_opt,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(residual_function, parameters_init, [], [], options); toc; -percent = (param_opt*100./[getParams_Cabbage();2.48252;0;0])-100 +percent = (param_opt*100./[getParamsCabbage();2.48252;0;0])-100*param_opt param_opt -%% Joint confidence intervals -%CI = nlparci(param_opt,residual,'jacobian',jacobian); +%% Joint confidence intervals +% if Statistics&ML toolbox is installed, use nlparci instead of helper: +% CI = nlparci(param_opt,residual,'jacobian',jacobian); +CI = computeCI(param_opt, residual, jacobian, 'alpha', 0.05); %% Solution with estimated parameters initialvalues_opt = param_opt(10:12); diff --git a/toolbox/examples/predatorpreyFilippov/runPredatorPrey.m b/toolbox/examples/predatorpreyFilippov/runPredatorPrey.m index b8974100..eea2a0ee 100644 --- a/toolbox/examples/predatorpreyFilippov/runPredatorPrey.m +++ b/toolbox/examples/predatorpreyFilippov/runPredatorPrey.m @@ -5,7 +5,7 @@ % Integrator intIfdiff = @ode45; intOptions = odeset('reltol', 1e-5, 'abstol', 1e-12); -eulerStep = 1e-7; +eulerStep = 1e-7; % Parameter values required for Shilnikov behavior in the paper (see RHS file). m = 0.790; diff --git a/toolbox/examples/sign/sign_inconsistent_rhs.m b/toolbox/examples/sign/rhsSignInconsistent.m similarity index 76% rename from toolbox/examples/sign/sign_inconsistent_rhs.m rename to toolbox/examples/sign/rhsSignInconsistent.m index 47b3c402..42a78343 100644 --- a/toolbox/examples/sign/sign_inconsistent_rhs.m +++ b/toolbox/examples/sign/rhsSignInconsistent.m @@ -1,4 +1,4 @@ -function dx = sign_inconsistent_rhs(t,x,p) +function dx = rhsSignInconsistent(t,x,p) dx = zeros(1,1); dx(1) = 2 - 4*sign(x(1)-p(1)); diff --git a/toolbox/examples/sign/sign_inconsistent_rhs_2d.m b/toolbox/examples/sign/rhsSignInconsistent2d.m similarity index 78% rename from toolbox/examples/sign/sign_inconsistent_rhs_2d.m rename to toolbox/examples/sign/rhsSignInconsistent2d.m index be5ae881..39911eb2 100644 --- a/toolbox/examples/sign/sign_inconsistent_rhs_2d.m +++ b/toolbox/examples/sign/rhsSignInconsistent2d.m @@ -1,4 +1,4 @@ -function dx = sign_inconsistent_rhs_2d(t,x,p) +function dx = rhsSignInconsistent2d(t,x,p) dx = zeros(2,1); sign_value = sign(x(2)-p(1)); diff --git a/toolbox/examples/sign/sign_inconsistent_test.m b/toolbox/examples/sign/runSignInconsistent.m similarity index 92% rename from toolbox/examples/sign/sign_inconsistent_test.m rename to toolbox/examples/sign/runSignInconsistent.m index 62a067e9..3ccd5942 100644 --- a/toolbox/examples/sign/sign_inconsistent_test.m +++ b/toolbox/examples/sign/runSignInconsistent.m @@ -1,5 +1,5 @@ integrator = @ode45; -initstates = [1]; +initstates = 1; p = 0; doIfdiff = true; @@ -12,9 +12,9 @@ fprintf('Preprocessing...\n '); odeoptions = odeset( 'AbsTol', 1e-14, 'MaxStep', 2); - filename = 'sign_inconsistent_rhs'; + filename = 'rhsSignInconsistent'; datahandle = prepareDatahandleForIntegration(filename, ... - 'integrator', func2str(integrator), 'options', odeoptions); + 'integrator', integrator, 'options', odeoptions); fprintf('Integration with ifdiff/%s...\n ', func2str(integrator)); tic diff --git a/toolbox/examples/sign/sign_inconsistent_2d_test.m b/toolbox/examples/sign/runSignInconsistent2d.m similarity index 91% rename from toolbox/examples/sign/sign_inconsistent_2d_test.m rename to toolbox/examples/sign/runSignInconsistent2d.m index 52f7c72d..bca0ceff 100644 --- a/toolbox/examples/sign/sign_inconsistent_2d_test.m +++ b/toolbox/examples/sign/runSignInconsistent2d.m @@ -1,17 +1,16 @@ integrator = @ode45; t0 = 0; -tf = 1; -% tf = 0.5+30*eps; +tf = 1; % tf = 0.5+30*eps; timeinterval = [t0,tf]; initstates = [1, 1]; p = 0; fprintf('Preprocessing...\n '); odeoptions = odeset( 'AbsTol', 1e-14, 'MaxStep', 2); -filename = 'sign_inconsistent_rhs_2d'; +filename = 'rhsSignInconsistent2d'; tic datahandle = prepareDatahandleForIntegration(filename, ... - 'solver', func2str(integrator), 'options', odeoptions); + 'solver', integrator, 'options', odeoptions); toc fprintf('Integration with ifdiff/%s...\n ', func2str(integrator)); @@ -20,8 +19,6 @@ % sol_ode45 = ode45(@(t,x) sign_inconsistent_rhs_2d(t,x,p), timeinterval, initstates); toc - - % Read out number of switches during integration n_switches = max(size(datahandle.getData().SWP_detection.switchingpoints)); fprintf('Total number of switches during integration: %d \n', n_switches); diff --git a/toolbox/examples/spiral/rhsSpiral.m b/toolbox/examples/spiral/rhsSpiral.m new file mode 100644 index 00000000..953b37ea --- /dev/null +++ b/toolbox/examples/spiral/rhsSpiral.m @@ -0,0 +1,13 @@ +function dx = rhsSpiral(~, x, p) + % ODE showing a spiral behavior, while every cross of an axis results + % in a switch. In this way, the switching frequency goes to infinity. + dx = zeros(2,1); + threshold = 1e-3; + + if sign(x(1)) > threshold + dx(1) = -sign(x(1)) + 2*sign(x(2)); + dx(2) = -2*sign(x(1)) - sign(x(2)); + else + dx = [sign(x(1));sign(x(2))]; + end +end \ No newline at end of file diff --git a/toolbox/examples/spiral/spiral_test.m b/toolbox/examples/spiral/runSpiral.m similarity index 87% rename from toolbox/examples/spiral/spiral_test.m rename to toolbox/examples/spiral/runSpiral.m index 1be3e5e8..19e27571 100644 --- a/toolbox/examples/spiral/spiral_test.m +++ b/toolbox/examples/spiral/runSpiral.m @@ -6,10 +6,10 @@ p = []; fprintf('Preprocessing...\n '); -odeoptions = odeset( 'AbsTol', 1e-8, 'RelTol', 1e-6, 'MaxStep', 10000); -filename = 'spiral_rhs'; +odeoptions = odeset( 'AbsTol', 1e-8, 'RelTol', 1e-6); +filename = 'rhsSpiral'; tic -handlerhs_test = prepareDatahandleForIntegration(filename, 'integrator', func2str(integrator), 'options', odeoptions); +handlerhs_test = prepareDatahandleForIntegration(filename, 'integrator', integrator, 'options', odeoptions); toc % solution diff --git a/toolbox/examples/spiral/spiral_rhs.m b/toolbox/examples/spiral/spiral_rhs.m deleted file mode 100644 index 0b76ad8d..00000000 --- a/toolbox/examples/spiral/spiral_rhs.m +++ /dev/null @@ -1,8 +0,0 @@ -function dx = spiral_rhs(t, x, p) - % ODE showing a spiral behavior, while every cross of an axis results - % in a switch. In this way, the switching frequency goes to infinity. - dx = zeros(2,1); - - dx(1) = -sign(x(1)) + 2*sign(x(2)); - dx(2) = -2*sign(x(1)) - sign(x(2)); -end \ No newline at end of file diff --git a/toolbox/examples/statejump/bounceball/bounceballRHS.m b/toolbox/examples/statejump/bounceball/rhsBounceball.m similarity index 98% rename from toolbox/examples/statejump/bounceball/bounceballRHS.m rename to toolbox/examples/statejump/bounceball/rhsBounceball.m index a2864231..1337aad6 100644 --- a/toolbox/examples/statejump/bounceball/bounceballRHS.m +++ b/toolbox/examples/statejump/bounceball/rhsBounceball.m @@ -1,4 +1,4 @@ -function dx = bounceballRHS(~, x, p) +function dx = rhsBounceball(~, x, p) % basic bouncing ball. It loses energy on each bounce, but not by air friction. There is no % solution to Zeno's phenomenon when the ball is very low, just don't integrate that far. % p = [g, gamma] diff --git a/toolbox/examples/statejump/bounceball/bounceballZenoRHS.m b/toolbox/examples/statejump/bounceball/rhsBounceballZeno.m similarity index 96% rename from toolbox/examples/statejump/bounceball/bounceballZenoRHS.m rename to toolbox/examples/statejump/bounceball/rhsBounceballZeno.m index 06f7da8a..4a238355 100644 --- a/toolbox/examples/statejump/bounceball/bounceballZenoRHS.m +++ b/toolbox/examples/statejump/bounceball/rhsBounceballZeno.m @@ -1,4 +1,4 @@ -function dx = bounceballZenoRHS(~, x, p) +function dx = rhsBounceballZeno(~, x, p) % basic bouncing ball with Zeno's phenomenon treatment. It loses energy on each bounce, but % not by air friction. % p = [g, gamma, eByMassMin] diff --git a/toolbox/examples/statejump/bounceballError.m b/toolbox/examples/statejump/bounceballError.m index fd6f95eb..7d7dceae 100644 --- a/toolbox/examples/statejump/bounceballError.m +++ b/toolbox/examples/statejump/bounceballError.m @@ -94,7 +94,7 @@ standardTic = tic; % Standard solution and sensitivities -datahandle = prepareDatahandleForIntegration('bounceballRHS', ... +datahandle = prepareDatahandleForIntegration('rhsBounceball', ... 'solver', func2str(solverStandard), ... 'options', optionsStandard); solStandard = solveODE(datahandle, [t0 tEnd], x0, p); @@ -130,7 +130,7 @@ % Strict solution and sensitivities strictTic = tic; % Strict solution and sensitivities -datahandle = prepareDatahandleForIntegration('bounceballRHS', ... +datahandle = prepareDatahandleForIntegration('rhsBounceball', ... 'solver', func2str(solverStrict), ... 'options', optionsStrict); solStrict = solveODE(datahandle, [t0 tEnd], x0, p); diff --git a/toolbox/examples/statejump/bounceballZeno.m b/toolbox/examples/statejump/bounceballZeno.m index 8b9b8d99..7fd4b66b 100644 --- a/toolbox/examples/statejump/bounceballZeno.m +++ b/toolbox/examples/statejump/bounceballZeno.m @@ -14,8 +14,8 @@ tEnd = 25; % zeno begins (analytically) at 20.3935964107270316 for these particular data x0 = [h0; v0]; -datahandle = prepareDatahandleForIntegration('bounceballZenoRHS', ... - 'solver', func2str(solver), ... +datahandle = prepareDatahandleForIntegration('rhsBounceballZeno', ... + 'solver', solver, ... 'options', options); sol = solveODE(datahandle, [t0 tEnd], x0, p); diff --git a/toolbox/examples/statejump/jumpSensitivities.m b/toolbox/examples/statejump/jumpSensitivities.m index 32b2507b..79b8bbed 100644 --- a/toolbox/examples/statejump/jumpSensitivities.m +++ b/toolbox/examples/statejump/jumpSensitivities.m @@ -42,8 +42,8 @@ % IFDIFF solution (standard tolerances) datahandle = prepareDatahandleForIntegration( ... - 'jumpSensitivityRHS', ... - 'solver', func2str(solverStandard), ... + 'rhsJumpSensitivity', ... + 'solver', solverStandard, ... 'options', optionsStandard); solStandard = solveODE(datahandle, [t0 tEnd], x0, p); @@ -68,8 +68,8 @@ % IFDIFF solution with strict tolerances. The smaller FDstep makes most of the difference, % but the strict ODE tolerances also help a little bit datahandle = prepareDatahandleForIntegration( ... - 'jumpSensitivityRHS', ... - 'solver', func2str(solverStrict), ... + 'rhsJumpSensitivity', ... + 'solver', solverStrict, ... 'options', optionsStrict); solStrict = solveODE(datahandle, [t0 tEnd], x0, p); diff --git a/toolbox/examples/statejump/jumpSensitivities/jumpSensitivityRHS.m b/toolbox/examples/statejump/jumpSensitivities/rhsJumpSensitivity.m similarity index 89% rename from toolbox/examples/statejump/jumpSensitivities/jumpSensitivityRHS.m rename to toolbox/examples/statejump/jumpSensitivities/rhsJumpSensitivity.m index 6d121df6..8db37471 100644 --- a/toolbox/examples/statejump/jumpSensitivities/jumpSensitivityRHS.m +++ b/toolbox/examples/statejump/jumpSensitivities/rhsJumpSensitivity.m @@ -1,4 +1,4 @@ -function dx = jumpSensitivityRHS(~, x, p) +function dx = rhsJumpSensitivity(~, x, p) %JUMPSENSITIVITYRHS A simple jumping model for testing sensitivity computation across jumps. dx = 0; sigma = x(1) - (1/p(1)); diff --git a/toolbox/examples/statejump/bounceball.m b/toolbox/examples/statejump/runBounceball.m similarity index 90% rename from toolbox/examples/statejump/bounceball.m rename to toolbox/examples/statejump/runBounceball.m index f0295246..be415ef8 100644 --- a/toolbox/examples/statejump/bounceball.m +++ b/toolbox/examples/statejump/runBounceball.m @@ -12,8 +12,8 @@ tEnd = 20; % zeno begins (analytically) at 20.3935964107270316 for these particular data x0 = [h0; v0]; -datahandle = prepareDatahandleForIntegration('bounceballRHS', ... - 'solver', func2str(solver), ... +datahandle = prepareDatahandleForIntegration('rhsBounceball', ... + 'solver', solver, ... 'options', options); sol = solveODE(datahandle, [t0 tEnd], x0, p); diff --git a/toolbox/examples/statejump/sinCosMain.m b/toolbox/examples/statejump/sinCosMain.m index 1fd53d1c..d6fe83f6 100644 --- a/toolbox/examples/statejump/sinCosMain.m +++ b/toolbox/examples/statejump/sinCosMain.m @@ -6,7 +6,7 @@ p = 0; x0 = [0; 1]; -datahandle = prepareDatahandleForIntegration('sinCosRHS', 'solver', func2str(solver), 'options', options); +datahandle = prepareDatahandleForIntegration('sinCosRHS', 'solver', solver, 'options', options); sol = solveODE(datahandle, [t0 tF], x0, p); T = linspace(t0, tF, 100); diff --git a/toolbox/examples/subway/Sensitivities/Plot_Sensitivities_NewYork.m b/toolbox/examples/subway/Sensitivities/Plot_Sensitivities_NewYork.m index 0eeb8ff5..9955fe06 100644 --- a/toolbox/examples/subway/Sensitivities/Plot_Sensitivities_NewYork.m +++ b/toolbox/examples/subway/Sensitivities/Plot_Sensitivities_NewYork.m @@ -5,7 +5,7 @@ parameters = nysscc_getPhysicsParameters_vector(); odeoptionssubwaymodel = odeset( 'AbsTol', 1e-20, 'RelTol', 1e-10); -datahandle = prepareDatahandleForIntegration('newYorkCitySubwayModelRhs_wrapped', 'integrator', func2str(integrator), 'options', odeoptionssubwaymodel); +datahandle = prepareDatahandleForIntegration('rhsNewYorkCitySubwayModelWrapped', 'integrator', integrator, 'options', odeoptionssubwaymodel); sol = solveODE(datahandle, tspan, initValues, parameters); sol_ode45 = ode45(@(t,x) newYorkCitySubwayModelRhs(t,x, nysscc_getPhysicsParameters()),tspan, initValues,odeoptionssubwaymodel); diff --git a/toolbox/examples/subway/Sensitivities/newYorkCitySubwayModelRhs_wrapped.m b/toolbox/examples/subway/Sensitivities/rhsNewYorkCitySubwayModelWrapped.m similarity index 98% rename from toolbox/examples/subway/Sensitivities/newYorkCitySubwayModelRhs_wrapped.m rename to toolbox/examples/subway/Sensitivities/rhsNewYorkCitySubwayModelWrapped.m index 6e6fd0a1..da396b08 100644 --- a/toolbox/examples/subway/Sensitivities/newYorkCitySubwayModelRhs_wrapped.m +++ b/toolbox/examples/subway/Sensitivities/rhsNewYorkCitySubwayModelWrapped.m @@ -1,4 +1,4 @@ -function dy = newYorkCitySubwayModelRhs_wrapped(t, y, p) +function dy = rhsNewYorkCitySubwayModelWrapped(t, y, p) pp.Tmax = p(1); pp.S = p(2); diff --git a/toolbox/examples/subway/newYorkCitySubwayModelRhs.m b/toolbox/examples/subway/rhsNewYorkCitySubwayModel.m similarity index 98% rename from toolbox/examples/subway/newYorkCitySubwayModelRhs.m rename to toolbox/examples/subway/rhsNewYorkCitySubwayModel.m index e51ec295..8b96e548 100644 --- a/toolbox/examples/subway/newYorkCitySubwayModelRhs.m +++ b/toolbox/examples/subway/rhsNewYorkCitySubwayModel.m @@ -1,7 +1,7 @@ % newYorkCitySubwayRhs02.m -> version two % right-hand side function of the new york city subway model -function dy = newYorkCitySubwayModelRhs(t, y, p) +function dy = rhsNewYorkCitySubwayModel(t, y, p) % 't': time % 'y': states % 'p': parameters diff --git a/toolbox/examples/subway/newYorkCityExampleIntegrationForThesis.m b/toolbox/examples/subway/runNewYorkCityExampleIntegrationForThesis.m similarity index 97% rename from toolbox/examples/subway/newYorkCityExampleIntegrationForThesis.m rename to toolbox/examples/subway/runNewYorkCityExampleIntegrationForThesis.m index b3f8fa0e..d7ca2a72 100644 --- a/toolbox/examples/subway/newYorkCityExampleIntegrationForThesis.m +++ b/toolbox/examples/subway/runNewYorkCityExampleIntegrationForThesis.m @@ -1,15 +1,14 @@ % new York city integration integrator = @ode45; - timeinterval = [0.0, 65.0]; initstates = [0.0, 0.0, 0.0].'; odeoptionssubwaymodel = odeset( 'AbsTol', 1e-20, 'RelTol', 1e-10); %% tic -hdlnyc = prepareDatahandleForIntegration('newYorkCitySubwayModelRhs', ... - 'integrator', func2str(integrator), 'options', odeoptionssubwaymodel); +hdlnyc = prepareDatahandleForIntegration('rhsNewYorkCitySubwayModel', ... + 'integrator', integrator, 'options', odeoptionssubwaymodel); toc %% @@ -22,7 +21,7 @@ toc -sol = ode45(@(t,x) newYorkCitySubwayModelRhs(t,x, nysscc_getPhysicsParameters()),timeinterval, initstates,odeoptionssubwaymodel); +sol = ode45(@(t,x) rhsNewYorkCitySubwayModel(t,x, nysscc_getPhysicsParameters()),timeinterval, initstates,odeoptionssubwaymodel); %% hold on