MATLAB Answers

0

How to declare time dependent symbolic array?

Asked by Marlon Saveri Silva on 31 Mar 2018
Latest activity Edited by Walter Roberson
on 14 Oct 2018 at 1:03

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.

  0 Comments

Sign in to comment.

3 Answers

Answer by Marlon Saveri Silva on 31 Mar 2018

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

  4 Comments

Show 1 older comment

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

 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

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

Sign in to comment.


Answer by Sven
on 4 Oct 2018 at 17:10
Edited by Sven
on 4 Oct 2018 at 19:38

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

  0 Comments

Sign in to comment.


Answer by ahmad mahmoodi on 13 Oct 2018 at 18:23
Edited by Walter Roberson
on 13 Oct 2018 at 18:59

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!

  2 Comments

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.

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

Sign in to comment.