diff --git a/+misc/solveGenericBvpFirstOrder.m b/+misc/solveGenericBvpFirstOrder.m index 44c7ef42b1996a31f37ade431b5015bfbd359f53..9a3e867979b3ded72881e23b0e91eda382e5dd4d 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