Why is my variable being changed without assignment

Hello,

I have searched quite a lot on Google and cannot find anything related to what is happening, and just in general programming terms I cannot understand it either. I'm stumped and hope someone can shed some light. An analogy of what seems to be happening is this:

A = 1

B = A + 2 = 3 (and now A = 3!!!)

This makes no sense to me, how can the original value of A be getting changed without. This is my first time using Live Script and symbols, so perhaps there is something different around those concepts.

My code is below. The issue is that when I substitute 0 into yh(zero), it is also changing the value of yh. Thus when I perform a subs to calc yhNat (last line), the value it is referencing is no longer yh but actually yh(zero). If I leave out the calc of yh(zero) then yhNat calculates perfectly (15*exp((-2*t)) + 16*exp((-3*t))). The only thing I can think of is because the initial value assigned to both yh and yh(zero) is the same, perhaps Matlab is conflating/linking them somehow. Originally I just made yh(zero) a subs of yh, replacing all t's with 0, but that didn't work, hence the workaround.

Any help would be greatly appreciated, as I am fully stuck and all my years of messing around with programming languages aren't helping at all to understand why this is happening.

% Homogenous solution
InitialVoltage = 1.5;
InitialCurrent = 2;
c_UD = 1/6;
%from earlier code roots.1 = -3, roots.2 = -2

syms yh yh(zero) yhNat dyh(t) dyh(zero) c_1 c_2

yh = c_1*exp(roots(1)*t) + c_2*exp(roots(2)*t)
dyh(t) = diff(yh)
dyh(zero) = subs(dyh(t),[t],[0])  == InitialCurrent/c_UD

%Evaluate yh(t) at t = 0:
yh(zero) = c_1*exp(roots(1)*t) + c_2*exp(roots(2)*t);      %######### These two lines break it
yh(zero) = subs(yh(zero), [t], [0]) == vpa(InitialVoltage) %######### These two lines break it

%coefficients = solve([dyh(zero), yh(zero)], [c_1, c_2])

%Natural response
yhNat = subs(yh ,[c_1, c_2], [16, 15])

[Paul]

Hi Robyn

I'm going to guess that yh(zero) isn't what you want it to be. Consider

syms yh yh(zero)
whos y* z*

  Name      Size            Bytes  Class     Attributes

  yh        1x1                 8  symfun              
  zero      1x1                 8  sym                 

At this point, yh is an unassigned symfun object with sym object named zero as its argument. Based on the code, it looks like you were expectin yh and yh(zero) to be different things, but as we can see that's not the case.

Also, I'm not sure what that expression involving roots is supposed to be. As we can see

roots(1)

ans =

  0x1 empty double column vector

is probably not what's needed.

It will be easier for someone to help of you show the mathematical equations that you're tryting to implement. You can upload a screenshot using the Paperclip icon on the Insert menu. Of you can use the Sigma icon to edit your post with formatted equations if you're proficient with latex.

[Robyn Galliers]

Hi Paul,

Thanks for your reply. Yes, I very much expected yh and yh(zero) to be different variables.I thoughts symbols were a bunch of characters that were treated as one entity, i.e. brackets wouldn't be counted.

I am trying to recreate a worked example from my textbook to help me understand what is going on and then later build upon. It is the natural response of a RLC circuit, i.e. a 2nd order differential equation. An excerpt from the book is attached.

As for the roots bit, I am using an equation solver earlier in my live script to calculate them but didn't want to include all my code. But essentially it is doing

syms s;
CharEqn = s^2 + (r_UD/l_UD)*s + (1/(l_UD*c_UD)) == 0;
roots = solve(CharEqn,s,"MaxDegree",4);

[Robyn Galliers]

Paul,

I just changed my variable from yh(zero) to yhzero, and it all suddenly works, so you were bang on the money. The hours I wasted yesterday pleading with my computer to make sense to me, and it was all because of that :(

Thank you for your help, I greatly appreciate it!

Answer 1

[Sam Chak]

Hi @Robyn Galliers

Considering the given 1st-order differential equation with the parameters solved by @Sulaymon Eshkabilov, are the following responses expected? If that's the case, then it appears that this is indeed a 2nd-order system since it possesses two eigenvalues (a technical term for the roots in this context).

root    = [-3, -2];
c1      = -15.0;
c2      =  16.5;
dyh     = @(t, y) root(1)*c1*exp(root(1)*t) + root(2)*c2*exp(root(2)*t);
t       = 0:0.01:6;
InitVol = 1.5;
c_UD    = 1/6;
InitCur = dyh(0, InitVol)*c_UD

InitCur = 2

sol     = ode45(dyh, t, InitVol);
[y, yp] = deval(sol, t);

figure(1)
subplot(211)
plot(t, y),  grid on, xlabel('Time'), ylabel('Voltage')
subplot(212)
plot(t, yp*c_UD), grid on, xlabel('Time'), ylabel('Current')

Comparison with the simulation of a 2nd-order system:

%% State-space system
A   = [0, 1; -(root(1)*root(2)), (root(1) + root(2))];
B   = [0; 1];
C   = diag([1, c_UD]);
D   =  0;
sys = ss(A, B, C, D)

sys =
 
  A = 
       x1  x2
   x1   0   1
   x2  -6  -5
 
  B = 
       u1
   x1   0
   x2   1
 
  C = 
           x1      x2
   y1       1       0
   y2       0  0.1667
 
  D = 
       u1
   y1   0
   y2   0
 
Continuous-time state-space model.

ev  = eig(sys)      % check if the system has the same eigenvalues as the roots

ev = 2x1

-2.0000 -3.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

x0  = [InitVol; InitCur/c_UD];
u   = zeros(1, numel(t));
figure(2)
lsim(sys, u, t, x0), grid on

[Robyn Galliers]

Hi Sam,

Wow, thanks for your reply. You're using code I don't recognise, but you have drawn the graph I was hoping to get to. Sulyamon's response doesn't fix my issue, however Paul's comment highlighted that the symbols I am using aren't distinct like I thought they were. Becuase what I am trying to write is a step-by- step worked example, I am going to stick with that approach, but thank you for your time in replying and showing a different way of getting to the graph I was after.

[Sam Chak]

Hi @Robyn Galliers

I understand what you're looking for. If you want to display the analytical solution in a format similar to what is shown in the textbook, you can utilize the 'dsolve()' command, which is the standard method in MATLAB for symbolically finding the analytical solution of a differential equation.

However, the approach demonstrated by @Sulaymon Eshkabilov follows the traditional pen-and-paper procedure for solving the differential equation, just like in the textbook.

syms y(t)

%% parameters
InitVol = 1.5;
c_UD    = 1/6;
InitCur = 2;

%% setup differential equation and initial values
Dy      = diff(y,t);
DiffEqn = diff(y,t,2) + 5*Dy + 6*y == 0;
y0      = [y(0)==InitVol, Dy(0)==InitCur/c_UD];

%% Call dsolve
ySol(t) = dsolve(DiffEqn, y0);
ySol(t) = expand(ySol);
ySol(t) = vpa(ySol)

ySol(t) = 

[Sam Chak]

@Robyn Galliers: You're using code I don't recognise...

No need to worry! The state-space method (utilizing ss, eig, lsim) is typically covered in Year 3 of the curriculum. However, "Signals and Systems" is usually taught in Year 2.

From the 2nd-order differential equation

, where the external input signal (on right-hand side) is 0,

it can be rearranged to become

.

Next, let the state variables be

,

then we have

.

.

Rewrite the equations in matrix form:

State equation

Output equation

Answer 2

[Sulaymon Eshkabilov]

In your provided code, a couple of inconsistencies.

A first one is that roots()) is a MATLAB's builtin function. And if you want to assign a variable root, then it should be done in this way:

root = [-3, -2];

% It can be proceed:
% Homogenous solution
InitialVoltage = 1.5;
InitialCurrent = 2;
c_UD = 1/6;
% Assign the variable root values: 
root= [-3, -2];

syms yh yh(zero) yhNat dyh(t) dyh(zero) c_1 c_2

yh = c_1*exp(root(1)*t) + c_2*exp(root(2)*t)

yh = 

dyh(t) = diff(yh, t)

dyh(t) = 

dyh(zero) = subs(dyh(t),t,0)  == InitialCurrent/c_UD

dyh(zero) = 

% Evaluate yh(t) at t = 0:
yh(zero) = c_1*exp(root(1)*t) + c_2*exp(root(2)*t);       %######### These two lines break it
yh(zero) = subs(yh(zero), t, 0) == vpa(InitialVoltage) %######### These two lines break it

yh(zero) = 

Note how the solutions are found and validated:

% Solutions:
coefficients = solve([dyh(zero), yh(zero)], [c_1, c_2])

coefficients = struct with fields:
    c_1: -15.0
    c_2: 16.5

Solution = [coefficients.c_1, coefficients.c_2]

Solution = 

% Validation:

yhNat = subs(yh ,[c_1, c_2], Solution)

yhNat(zero) = 

Math works :)

All the best!

[Robyn Galliers]

Hi there,

Thanks for your reply. I have just explained about the roots in my response to the previous reply, I am using a solver to calculate the roots. But good point on roots() being a built-in function. I also attached to that response a screenshot of the textbook example I am working through.

The problem with the code you have above is that the final line is stil equating yhNat(zero) to 1.5, when it should be 16.5exp(-2t) - 15exp(-3t). In my mind it should take yh (c1exp(-3t)+c2exp-2t) and subs in the calculated coefficients values. But for some reason it is taking yh(zero) and subsing in the coefficient values.

The previous respondant suggested it was because yh and yh(zero) are not two different variables, so I might just have to pick different variables names. It's annoying as I wanted to try write them as close to the textbook as possible, but I wasted hours yesterday trying to write y_h(t), y_h(0), dy__h(t)/dt etc. Ideally I would be able to write the formulas exactly like in the textbook, but Matlab still be able to solve them.