Commit 4b67d09e authored by Jakob Gabriel's avatar Jakob Gabriel
Browse files

misc.solveGenericBvpFirstOrder: catch the case that jordan form is not...

misc.solveGenericBvpFirstOrder: catch the case that jordan form is not recognized due to numerical uncertainty -> round()
parent 5a7e535e
...@@ -53,6 +53,11 @@ E1 = I(:, diag(L0)>0); ...@@ -53,6 +53,11 @@ E1 = I(:, diag(L0)>0);
E2 = I(:, diag(L0)<0); E2 = I(:, diag(L0)<0);
% get jordan chains % get jordan chains
[V, J, lengthJb] = misc.jordan_complete(S); [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); mu = diag(J);
numJb = numel(lengthJb); numJb = numel(lengthJb);
Linv = inv(L); Linv = inv(L);
...@@ -106,6 +111,9 @@ for it = 1 : numJb ...@@ -106,6 +111,9 @@ for it = 1 : numJb
end % for kt = 1 : lenghtJb(it) end % for kt = 1 : lenghtJb(it)
end % for it = 1 : numJb end % for it = 1 : numJb
X = horzcat(X_ik{:}) / V; 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 if calculateW
W = horzcat(W_ik{:}) / V; W = horzcat(W_ik{:}) / V;
else else
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment