solveGenericBvpFirstOrder.m

function [X, W] = solveGenericBvpFirstOrder(S, L, Az, A, A0, Q0, B, Cx, F, NameValue)
% MISC.SOLVEGENERICBVPFIRSTORDER solve a generic matrix boundary value problem (BVP)
% with 1st order derivatives.
%
% [X, W] = misc.solveGenericBvpFirstOrder(S, L, Az, A, A0, Q0, B, Cx, F, ...
%											"Cw", Cw, "C", C, "D", D, "E", E)
% returns the solution of the BVP
%	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)
% [X] = solveGenericBvpFirstOrder(S, L, Az, A, A0, Q0, B, Cx, F)
% returns the solution of the BVP
%	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
%
% The parameter L, A may be constant matrices or quantity.Discrete.
% Meanwhile, S, Q0, B, Cw, F, C, E are constant matrices and Cx is a model.Output 
% object. The optional parameter Cw, C, D, E must be given as a name-value-pair.

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

% Check if W has to be calculated or not
if isempty(fieldnames(NameValue))
	calculateW = false;
	Cw = [];
else
	calculateW = true;
	Cw = NameValue.Cw;
	C = NameValue.C;
	D = NameValue.D;
	E = NameValue.E;
end

L0 = L.at(0);
I = eye(size(L0));
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);

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

numerator = cell(1, numJb);
for it = 1 : numJb
	my = mu(sum(lengthJb(1:(it-1)))+1);
	Psi = expm(cumInt(...
			Linv.subs("z", "eta") * (I * my - Az.subs("z", "eta")), "eta", "zeta", "z"));
	%M = Psi(z, 0, my) * (E1 + E2 * Q0) ...
	%	- int_0^z Psi(z, zeta, my) inv(L(zeta)) * A0(zeta) dzeta
	M = Psi.subs("zeta", 0) * (E1 + E2 * Q0) ...
		- cumInt(Psi * subs(Linv * A0, "z", "zeta"), "zeta", 0, "z");
	numerator{it} = Cx.out(M);
	if calculateW
		numerator{it} = numerator{it} - (Cw / (my * eye(size(C)) - C)) * E;
	end
	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;
		if calculateW
			D_ik = D * v_ik;
		end
		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) ...
		%				( X_i(k-1)(zeta) + A_ik(zeta) )	dzeta
		m_ik = Psi.subs("zeta", 0) * E2 * B_ik ...
				+ cumInt(Psi * subs(Linv * (X_ikMinus1 + A_ik), "z", "zeta"), "zeta", 0, "z");
		
		E1TX_ik0 = numerator{it} \ (F_ik - Cx.out(m_ik));
		if calculateW
			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);
		end
		X_ik{columnOfV} = M * E1TX_ik0 + m_ik;
	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
	W = [];
end % if calculateW
end % function solveGenericBvpFirstOrder()