How to declare time dependent symbolic array?

Hello,

I need to solve the following differential equation

odes=diff(T) == A*T + B;
[VF,Subs] = odeToVectorField(odes);
Sys = matlabFunction(VF,'Vars',{'t','Y'});
[T,Y] = ode45(Sys, [0 5], T0);

Therefore, I need to declare T(t). It could be done as below

syms T1(t) T2(t) T3(t) T4(t) T5(t)

However, the the length of T vector is not necessarily 5, but an user defined number. So, I tried the following:

syms t; T(t)=sym('T',[Num,1]);

But T was not considered as a t function and. When "odes" differentiates it, dT/dt is considered zero and I can't solve the problem.

So, how to create time dependent symbolic variables in an array with Num lines?

This question complements another.

Answer 1

[Marlon Saveri Silva]

I got a solution from another topic, but I'm not sure if it's the best, since some warnings appears on screen:

T=sym([]);
 for i=1:N
       syms(sprintf('T%d(t)', i)) %declare each element in the array as a single symbolic function
       TT = symfun(sym(sprintf('T%d(t)', i)), t); %declare each element to a symbolic "handle"
       T = [T;TT]; %paste the symbolic "handle" into an array 
 end

[Walter Roberson]

You should use cell arrays instead of [] arrays for the vector of functions. If you use [] arrays then you cannot select individual elements out of the vector.

[Q]

Hi Walter, could you please give an example of how to use the cell arrays for the above problem? Thanks.

[Walter Roberson]

 T=cell(N,1);
 for I=1:N
       syms(sprintf('T%d(t)', i)) %declare each element in the array as a single symbolic function
       TT = symfun(sym(sprintf('T%d(t)', i)), t); %declare each element to a symbolic "handle"
       T{I} = TT; %paste the symbolic "handle" into an array 
 end

[Walter Roberson]

For R2018a and later,

 T = cell(N,1);
 for I=1:N
       TT = symfun(str2sym(sprintf('T%d(t)', i)), t); %declare each element to a symbolic "handle"
       T{I} = TT; %paste the symbolic "handle" into an array 
 end

Answer 2

[Sven]

Hi Marlon,

if you substitute sym with eval, defining T, then there is no warning anymore. Also, @Q please see below where I used a cell array.

for i=1:5
  syms(sprintf('T%d(t)', i))
  TT = symfun(eval(sprintf('T%d(t)', i)), t); % use eval instead of sym
  T{i} = TT; % use cell array
end

Answer 3

[ahmad mahmoodi]

Hello all,

I used the same method. However I can not multiply my cell array of coefficients (AB) into the cell array of symbolic functions (T, Tm). In the following is the abstract of my code:

T=cell(4*n,1);
Tm=cell(6*n,1);
AB=cell(4*n,6*n);
GF=cell(4*n,1);

for i=1:n

AB, Tm an GF (using your comments) and the relation between them are defined here

end

odes=diff(T)==AB*Tm+GF 
or
(odes=diff(T)==AB.*Tm+GF;)

errors:

Undefined operator '*' for input arguments of type 'cell'.
(Undefined operator '.*' for input arguments of type 'cell'.)

Also I tried to use matrix instead as follow:

T=sym(zeros(4*n,1));
Tm=sym(zeros(6*n,1));

for i=1:n

AB, Tm an GF using your comments as matrices and the relation between them are defined here

end

   odes=diff(T)==AB*Tm+GF;
   [VF,Subs] = odeToVectorField(odes);
   Sys = matlabFunction(VF,'Vars',{'t','Y'});
   Y0(1:4*n,1)=14.7;
   [T,Y] = ode45(Sys, [0 5], Y0);

and here is the error:

Error using mupadengine/feval (line 163)
Cannot convert the initial value problem to an equivalent dynamical system. Either the differential equations cannot be solved for the highest derivatives or
inappropriate initial conditions were specified.

Error in odeToVectorField>mupadOdeToVectorField (line 177)
T = feval(symengine,'symobj::odeToVectorField',sys,x,stringInput);

Error in odeToVectorField (line 126)
sol = mupadOdeToVectorField(varargin);

Error in Fb1 (line 154)
      [VF,Subs] = odeToVectorField(odes);

Thanks in advance for your reply!

[Walter Roberson]

It seems unlikely to me that Mathworks will ever define a meaning for the * operator between a cell array and anything else.

The point of using cell arrays was that the original poster was trying to define a vector of symbolic functions. MATLAB does not have any problem with vectors of symbolic expressions, but when you have vectors of symbolic functions, then the meaning of

TheVector(Value)

becomes ambiguous. Does it mean to index the vector at the given value, pulling out one of the functions? Or does it mean to execute all of the functions in turn, passing in the Value to each one? MATLAB decided that it should mean executing all of them in turn. Because of that, you cannot index into a vector of symbolic functions. A vector of symbolic functions becomes the same as if you had created a single symbolic function whose body was the vector.

When you are working with differential equations, chances are that what you need is a vector of symbolic expressions, rather than a vector of functions.

... Basically, you should probably be avoiding creating functions dynamically.

Cannot convert the initial value problem to an equivalent dynamical system. Either the differential equations cannot be solved for the highest derivatives or
inappropriate initial conditions were specified.

We will probably need your actual code (and the inputs) to analyze that.

[ahmad mahmoodi]

Hello Walter. Thank you for your quick response.

My problem consists of layers of differential equations. For the layers between each layer has 4 differential equations and 4 variables. The number of equations and the coefficients in the first and the last layer differ from the ones between (3 equations and 4 variable for the first layer and 4 equations 3 variables for the last layer). So totally I have n equations and n variables. Also two of the variables from each that layer are used in the next layer. So I have a for loop that I need to use if loop inside it. All of my variables are time dependant.

I solved the problem with just two layers by entering the equations manually using ode45 and the results match. But I need my code to calculate when I have multiple layers . Here is the code:

% parameters

c_p_grout=1500;                    
rho_grout=1460;                 
lambda_grout=1.5;               
c_p_fluid=4183;                 
rho_fluid=997.6;                
lambda_soil=2.51225;           
lambda_fluid=0.5913;            
m_dot=0.46;

L=100;                         
dl=50;                        
n=L/dl+1;                        
d_b=0.2;                        
d_a=0.04;                       
d_i=d_a-2*0.0037;               
d_s1=0.22;                      
d_z1=sqrt((d_b^2+d_s1^2)/2);    

R_conv=0.00670;                 
R_cond_1=0.08568;               
x=0.70275;                      
R_a=0.23701;                   
R_b=0.073;

R_g=4*R_b-R_conv-R_cond_1;     
R_gg=1/((1/(R_a-R_conv-R_cond_1-x*R_g))-(1/(1-x)/R_g))*4;
R_fg=R_conv+R_cond_1+x*R_g;
R_gb=(1-x)*R_g;

R_gg1=R_gg;
R_gg2=R_gg1;
R_fgred=R_fg/2;
R_gbred_wo_soil=R_gb/2;
R_ggred_wo_soil=1/((2/R_gg1)+(2/R_gg2));
R_bz1=log(d_z1/d_b)/(2*pi*lambda_soil);                           
R_gbred=R_gbred_wo_soil+2*R_bz1;
R_ggred=1/((1/R_ggred_wo_soil)+(1/2/R_gbred_wo_soil)-(1/2/R_gbred));
C_fluid=2*rho_fluid*c_p_fluid*pi/4*d_i^2*dl/2;  %C_fluid at layer 1
C_g=rho_grout*(pi/4)*((d_b^2/2)-(2*d_a^2))*c_p_grout*dl/2;            % In EES 2*d_a^2 instead of d_a^2
T_b=15; 
Tinlet=90; 

A=(-dl/2)*((1/C_g*R_ggred)+(1/C_g*R_gbred)+(1/C_g*R_fgred));
B=(dl/2)*(1/C_g*R_fgred);
C=(dl/2)*(1/C_g*R_ggred);
D=(-dl/2)*((1/C_fluid*R_fgred))-(m_dot*c_p_fluid/C_fluid);
E=(m_dot*c_p_fluid/C_fluid);
F=(dl/2)*(1/C_fluid*R_fgred);
G=(dl/2)*(T_b/C_g*R_gbred);
A1=2*A;
B1=2*B;
C1=2*C;
F1=2*F;
Gi=2*G;
D1=(-dl)*((1/C_fluid*R_fgred))-(m_dot*c_p_fluid/C_fluid);

T=sym(zeros(4*n,1));
Tm=sym(zeros(6*n,1));

for i=1:n

        syms(sprintf('T%d(t)', 4*i-3))
        TT=symfun(sym(sprintf('T%d(t)', 4*i-3)), t);
        T(4*i-3,:)=TT;
        syms(sprintf('T%d(t)', 4*i-2))
        TT1=symfun(sym(sprintf('T%d(t)', 4*i-2)), t);
        T(4*i-2,:)=TT1;
        syms(sprintf('T%d(t)', 4*i-1))
        TT2=symfun(sym(sprintf('T%d(t)', 4*i-1)), t);
        T(4*i-1,:)=TT2;                                     
        syms(sprintf('T%d(t)', 4*i))
        TT3=symfun(sym(sprintf('T%d(t)', 4*i)), t);
        T(4*i,:)=TT3;

        Tm(6*i-5,:)=TT;
        Tm(6*i-4,:)=TT1;
  %     Tm(6*i-3,:)    inside if loop
        Tm(6*i-2,:)=TT2;
        Tm(6*i-1,:)=TT3;
  %     Tm(6*i,:)      inside if loop

        if i==1

        %Tg1(i)=A*Tg1(i)+B*T1(i)+C*Tg2(i)+G;
        %Tg2(i)=A*Tg2(i)+B*T2(i)+C*Tg1(i)+G;
        %dT2(i)=D*T2(i)+E*T2(i+1)+F*Tg2(i);

        T(4*i-1,:)=Tinlet;
        Tm(6*i-3,:)=0;
        syms(sprintf('T%d(t)', 4*i+4))
        TT6=symfun(sym(sprintf('T%d(t)', 4*i+4)), t);
        Tm(6*i,:)=TT6;
        AB1=[A C 0 B 0 0; C A 0 0 B 0 ; 0 0 0 0 0 0; 0 F 0 0 D E];
        G1=[G; G; 0; 0];
        AB(4*(i-1)+1:4*(i-1)+4,6*(i-1)+1:6*(i-1)+6)=AB1;
        GF(4*(i-1)+1:4*(i-1)+4,1)=G1;

        elseif i==n

%     Tg1(i)=A*Tg1(i)+B*T1(i)+C*Tg2(i)+G;
%     Tg2(i)=A*Tg2(i)+B*T2(i)+C*Tg1(i)+G;
%     T1(i)=D*T1(i)+E2*T1(i-1)+F*Tg1(i);

        Tm(6*i,:)=0;
        syms(sprintf('T%d(t)', 4*i-5))
        TT7=symfun(sym(sprintf('T%d(t)', 4*i-5)), t);
        Tm(6*i-3,:)=TT7;
        Tm(6*i-1,:)=Tm(6*i-2,:);
        T(4*i-1,:)=T(4*i,:);
        AB3=[A C 0 B 0 0; C A 0 0 B 0 ; F 0 E D 0 0; 0 0 0 0 0 0];                      
        G3=[G; G; 0; 0];
        AB(4*(i-1)+1:4*(i-1)+4,6*(i-1)+1:6*(i-1)+6)=AB3;
        GF(4*(i-1)+1:4*(i-1)+4,1)=G3;

        else

%     Tg1(i)=A*Tg1(i)+B*T1(i)+C*Tg2(i)+G;
%     Tg2(i)=A*Tg2(i)+B*T2(i)+C*Tg1(i)+G;
%     T1(i)=D*T1(i)+E2*T1(i-1)+F*Tg1(i);
%     T2(i)=D*T2(i)+E2*T2(i+1)+F*Tg2(i)

        syms(sprintf('T%d(t)', 4*i-5))
        TT8=symfun(sym(sprintf('T%d(t)', 4*i-5)), t);
        Tm(6*i-3,:)=TT8;
        syms(sprintf('T%d(t)', 4*i+4))
        TT9=symfun(sym(sprintf('T%d(t)', 4*i+4)), t);
        Tm(6*i,:)=TT9;
        AB2=[A1 C1 0 B1 0 0; C1 A1 0 0 B1 0 ; F1 0 E D1 0 0; 0 F1 0 0 D E]; 
        G2=[Gi; Gi; 0; 0];
        AB(4*(i-1)+1:4*(i-1)+4,6*(i-1)+1:6*(i-1)+6)=AB2;
        GF(4*(i-1)+1:4*(i-1)+4,1)=G2;

        end 

end

        odes=diff(T)==AB*Tm+GF;
        [VF,Subs] = odeToVectorField(odes);
        Sys = matlabFunction(VF,'Vars',{'t','Y'});
        Y0(1:4*n,1)=14.7;
        [T,Y] = ode45(Sys, [0 5], Y0);

I also tried to solve using

function dT = SRBZ(t, T)

the previous loops and  coeficients here... 

end

T0(1:4*n,1)=14.7;
[t,T] = ode45(@SRBZ,[0 2],T0);

But it gives error again and can not solve

Answer 4

[Andrés Arciniegas]

Greetings. This simple loop generates a vector of time dependent variables. Works smoothly and might be useful to someone :).

Addressing the problem above, the script generates the vector: .

% Script to create time dependent variables
Tt =[] % Array
syms(var_ind)
for i=1:n
   var_str = strcat('T',num2str(i),'(t)'); % Builds the symbolic variable string
   syms(var_str); 
   Tt = [Tt str2sym(var_str)]; % store the variable in the array
end