From 4b67d09ec5e612b820fc722aecb3eb2eef6ccffc Mon Sep 17 00:00:00 2001
From: Jakob Gabriel <jakob.gabriel@uni-ulm.de>
Date: Tue, 20 Oct 2020 15:02:20 +0200
Subject: [PATCH] misc.solveGenericBvpFirstOrder: catch the case that jordan
form is not recognized due to numerical uncertainty -> round()
---
+misc/solveGenericBvpFirstOrder.m | 8 ++++++++
1 file changed, 8 insertions(+)
diff --git a/+misc/solveGenericBvpFirstOrder.m b/+misc/solveGenericBvpFirstOrder.m
index 44c7ef4..9a3e867 100644
--- a/+misc/solveGenericBvpFirstOrder.m
+++ b/+misc/solveGenericBvpFirstOrder.m
@@ -53,6 +53,11 @@ E1 = I(:, diag(L0)>0);
E2 = I(:, diag(L0)<0);
% get jordan chains
[V, J, lengthJb] = misc.jordan_complete(S);
+if abs(det(V)) < 1e-12
+ % in some badly scaled cases, the jordan blocks of S are not calculated correctly, which leads
+ % to badly scaled modal matrices V. Hence, better results are obtained after round().
+ [V, J, lengthJb] = misc.jordan_complete(round(S, 15));
+end
mu = diag(J);
numJb = numel(lengthJb);
Linv = inv(L);
@@ -106,6 +111,9 @@ for it = 1 : numJb
end % for kt = 1 : lenghtJb(it)
end % for it = 1 : numJb
X = horzcat(X_ik{:}) / V;
+assert(MAX(abs(imag(X))) < 1e-6, "there is a imaginary part with a max abs value of " ...
+ + num2str(MAX(abs(imag(X)))) + " in the solution");
+X = real(X);
if calculateW
W = horzcat(W_ik{:}) / V;
else
--
GitLab