solveGenericBvpFirstOrder.m 3.77 KB
Newer Older
1
function [X, W] = solveGenericBvpFirstOrder(S, L, Az, A, A0, Q0, B, Cx, F, NameValue)
2
3
4
% MISC.SOLVEGENERICBVPFIRSTORDER solve a generic matrix boundary value problem (BVP)
% with 1st order derivatives.
%
5
% [X, W] = misc.solveGenericBvpFirstOrder(S, L, Az, A, A0, Q0, B, Cx, F, ...
6
%											"Cw", Cw, "C", C, "D", D, "E", E)
7
% returns the solution of the BVP
8
9
10
11
%	L(z) X(z)' + Az(z) X(z) + A0(z) E1.' X(0) - X(z) S = A(z)
%								 (E2.' - Q0 E1.') X(0) = B
%										  Cx[X] + Cw W = F
%											 C W - W S = D + E E1.' X(0)
12
% [X] = solveGenericBvpFirstOrder(S, L, Az, A, A0, Q0, B, Cx, F)
13
% returns the solution of the BVP
14
15
16
%	L(z) X(z)' + Az(z) X(z) + A0(z) E1.' X(0) - X(z) S = A(z)
%								 (E2.' - Q0 E1.') X(0) = B
%												 Cx[X] = F
17
18
%
% The parameter L, A may be constant matrices or quantity.Discrete.
19
% Meanwhile, S, Q0, B, Cw, F, C, E are constant matrices and Cx is a model.Output 
20
21
% object. The optional parameter Cw, C, D, E must be given as a name-value-pair.

22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
arguments
	S (:, :) double;
	L (:, :) quantity.Discrete;
	Az (:, :) quantity.Discrete;
	A (:, :) quantity.Discrete;
	A0 (:, :) quantity.Discrete;
	Q0 (:, :) double;
	B (:, :) double;
	Cx (:, :) model.Output;
	F (:, :) double;
	NameValue.Cw (:, :) double;
	NameValue.C (:, :) double;
	NameValue.D (:, :) double;
	NameValue.E (:, :) double;
end
37

38
% Check if W has to be calculated or not
39
if isempty(fieldnames(NameValue))
40
	calculateW = false;
41
	Cw = [];
42
43
else
	calculateW = true;
44
45
46
47
	Cw = NameValue.Cw;
	C = NameValue.C;
	D = NameValue.D;
	E = NameValue.E;
48
49
end

50
L0 = L.at(0);
51
I = eye(size(L0));
52
53
54
55
E1 = I(:, diag(L0)>0);
E2 = I(:, diag(L0)<0);
% get jordan chains
[V, J, lengthJb] = misc.jordan_complete(S);
56
57
58
59
60
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
61
62
mu = diag(J);
numJb = numel(lengthJb);
63
Linv = inv(L);
64
65
66

X_ik = cell(1, numel(mu));
W_ik = cell(1, numel(mu));
67

68
69
70
numerator = cell(1, numJb);
for it = 1 : numJb
	my = mu(sum(lengthJb(1:(it-1)))+1);
71
	Psi = expm(cumInt(...
72
			Linv.subs("z", "eta") * (I * my - Az.subs("z", "eta")), "eta", "zeta", "z"));
73
	%M = Psi(z, 0, my) * (E1 + E2 * Q0) ...
74
	%	- int_0^z Psi(z, zeta, my) inv(L(zeta)) * A0(zeta) dzeta
75
	M = Psi.subs("zeta", 0) * (E1 + E2 * Q0) ...
76
		- cumInt(Psi * subs(Linv * A0, "z", "zeta"), "zeta", 0, "z");
77
78
	numerator{it} = Cx.out(M);
	if calculateW
79
		numerator{it} = numerator{it} - (Cw / (my * eye(size(C)) - C)) * E;
80
	end
81
82
83
84
85
	for kt = 1 : lengthJb(it)
		columnOfV = sum(lengthJb(1:(it-1)))+kt;
		v_ik = V(:, columnOfV);
		A_ik = A * v_ik;
		B_ik = B * v_ik;
86
87
88
		if calculateW
			D_ik = D * v_ik;
		end
89
90
91
92
93
94
95
96
97
98
99
		F_ik = F * v_ik;
		
		if kt > 1
			X_ikMinus1 = X_ik{columnOfV-1};
			W_ikMinus1 = W_ik{columnOfV-1};
		else
			X_ikMinus1 = zeros(size(L, 1), 1);
			W_ikMinus1 = zeros(size(Cw, 2), 1);
		end
		% m_ik = Psi(z, 0, my) E2 B_ik ...
		%		+ int_0^z Psi(z, zeta, my) L^-1(zeta) ...
100
		%				( X_i(k-1)(zeta) + A_ik(zeta) )	dzeta
101
		m_ik = Psi.subs("zeta", 0) * E2 * B_ik ...
102
				+ cumInt(Psi * subs(Linv * (X_ikMinus1 + A_ik), "z", "zeta"), "zeta", 0, "z");
103
		
104
		E1TX_ik0 = numerator{it} \ (F_ik - Cx.out(m_ik));
105
		if calculateW
106
107
108
			E1TX_ik0 = E1TX_ik0 + numerator{it} \ ...
					((Cw / (my * eye(size(C)) - C)) * (W_ikMinus1 + D_ik));
			W_ik{columnOfV} = (my * eye(size(C)) - C) \ (- E * E1TX_ik0 - W_ikMinus1 - D_ik);
109
		end
110
111
112
113
		X_ik{columnOfV} = M * E1TX_ik0 + m_ik;
	end % for kt = 1 : lenghtJb(it)
end  % for it = 1 : numJb
X = horzcat(X_ik{:}) / V;
114
115
116
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);
117
118
119
120
if calculateW
	W = horzcat(W_ik{:}) / V;
else
	W = [];
121
end % if calculateW
122
end % function solveGenericBvpFirstOrder()