Matrix differential equations, correct answer?

This is driving me crazy and I don't know another way to formulate the question but to type the code (It's really simple)

Step1: Create a matrix. (We will call it "D")

Step2: Solve a matrix differential equation.

Step3: Compute the result.

This is the code:

%%
clc;
clear all;
close all;

%% STEP 1. GET THE DATA
n_of_states = 28;

% Define the failure rates
dsw= 0.0215;
dd= 0.005186;
daux_gen=0.262+0.161+0.03+0.012;

dc= [0.1294 0.4334  0.8725];
dar= 6*dsw+daux_gen;
dpr= 6*dd;

%% STEP 2. CREATION OF A MATRIX (transition intensity matrix)
%The matrix dimension will be: (n_of_states x n_of_states)
D=zeros(n_of_states);

%We introduce the failure rate of the active rectifier
D(1,13)=dar;
D(4,14)=dar;
D(7,15)=dar;
D(10,16)=dar;
D(2,17)=dar;
D(5,18)=dar;
D(8,19)=dar;
D(11,20)=dar;
D(3,21)=dar;
D(6,22)=dar;
D(9,23)=dar;
D(12,24)=dar;

%We introduce the failure rate of the capacitor (1)
D(1,2)=dc(1);
D(4,5)=dc(1);
D(7,8)=dc(1);
D(10,11)=dc(1);

%We introduce the failure rate of the capacitor (2)
D(2,3)=dc(2);
D(5,6)=dc(2);
D(8,9)=dc(2);
D(11,12)=dc(2);

%We introduce the failure rate of the capacitor (3)
D(3,25)=dc(3);
D(6,26)=dc(3);
D(9,27)=dc(3);
D(12,28)=dc(3);

%We introduce the failure rate of the passive rectifier
D(1,4)=3*dpr;
D(2,5)=3*dpr;
D(3,6)=3*dpr;

D(4,7)=2*dpr;
D(5,8)=2*dpr;
D(6,9)=2*dpr;

D(7,10)=dpr;
D(8,11)=dpr;
D(9,12)=dpr;

%Completition
vector_sum=sum(D,2); %We create a vector with the sum of each row

for i = 1:n_of_states
    D(i,i)=-vector_sum(i);
end

Then (here it comes the problem), I do this:

syms p(t) [1 n_of_states];
Y=p;

%This is the matrix system
odes = diff(Y) == Y*D;

%This is the initial condition
C=zeros(1,n_of_states);
C(1)=1;
cond = p(0) == C;

%This should be the solution
[p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16,p17,p18,p19,p20,p21,p22,p23,p24,p25,p26,p27,p28] = dsolve(odes,cond);

reliability=p1+p2+p3+p4+p5+p6+p7+p8+p9+p10+p11+p12;
reliability=vpa(reliability,3);

and I obtain the function called "reliability".

Why is that function different to the one (called "rel") obtained doing this?

syms probabilidad(t) [1 n_of_states];
Y=probabilidad;

%This is the matrix system
ode = diff(Y) == Y*D;

%This is the initial condition
C=zeros(1,n_of_states);
C(1)=1;
condition = probabilidad(0) == C;

%This should be the solution
J = dsolve(ode,condition);
rel=J.probabilidad1+J.probabilidad2+J.probabilidad3+J.probabilidad4+J.probabilidad5+J.probabilidad6+J.probabilidad7+J.probabilidad8+J.probabilidad9+J.probabilidad10+J.probabilidad11+J.probabilidad12;
rel=vpa(rel,3);

And most importantly, WHICH ONE IS CORRECT??????? I'm desesperated!

Answer 1

[David Goodmanson]

Hi Sara,

this plot

figure(1)
tvec = 0:.01:6;
reliability_d = double(subs(reliability,t,tvec));
rel_d = double(subs(rel,t,tvec));
plot(tvec,reliability_d,tvec,rel_d)
grid on

shows that betwen the two versions, the first one has issues. As you can see, it's even increasing in the second half of the plot.

The second version 'rel' is correct, as verified here using an eigenvalue calculation. (In the last line of your code I upped the significant figures with

rel=vpa(rel,5)

). Due to the form of D, with all zeros in the lower left corner, time derivates of variables 1-12 do not depend on variables 13-28. Therefore the problem can be reduced to the upper left 12x12 corner of D.

tvec = 0:.01:6;
rel_d = double(subs(rel,t,tvec));         % your second result
% reduce problem to 12x12
D12 = D(1:12,1:12);
phi0 = C(1:12);                           % initial condition
[V lambda]= eig(D12);
phi0V = phi0*V;
inVSum = sum(inv(V),2);
cof = phi0V.*inVSum';                     % coefficients of exp(lambda*t) terms
[tvecx lambdax] = meshgrid(tvec,diag(lambda));  % time across,lambda down
phimat = exp(lambdax.*tvecx);
relnew = cof*phimat;

figure(2)
plot(tvec,relnew,tvec,rel_d +.01)        % add .01 for visual purposes
grid on

rel
[cof;diag(lambda)']                     % coefficients, lambdas

The last line shows that the solution has three exponentials of the form exp(lambda*t) with nonzero coefficients. If you scroll through the somewhat inconvenient rel line, you will find the same ones.

[Sara Linares]

Thank you very much!! You were really helpful!!