(ODE23) Unable to perform assignment because the left and right sides have a different number of elements error

Hello guys, I am coming across a problem when trying to formulate a differential equation to be simulated by ode23. My code iis as follows:

function dx=boost(t,x)
global v_in R C L T d
dx=zeros(2,1);
if pwm(t,T,d)==1
    dx(1) = v_in./L ;
    dx(2) = -x(2)./R*C ;
end
if pwm(t,T,d)==0 & x(1)>=0
    dx(1) = (v_in - x(2))./L ;
    dx(2) = x(1)./C - x(2)./R*C ;
end
if pwm(t,T,d)==0 & x(1)<0           %#ok<*AND2>
    dx(1) = 0 ;
    dx(2) = -x(2)./R*C ;
end
dx=[dx(1);dx(2)]; 
end

so the diff eq. takes three forms, as seen by the three if statements and they depend on the current outputs of another function. This other function is a pwm signal which controls a switch (this is being used to simulate a boost converter circuit but that is irrelevant).

the error 'Unable to perform assignment because the left and right sides have a different number of elements' refers to line 9 [ dx(1) = (v_in - x(2))./L ] . It makes no sense to me why I see this error since as far as I am aware both sides of the equations are scalar. Also I dont see this error for any other assigments in the function which all have similar forms.

Thank you.

[dpb]

Probably one of the globals is a vector...I notice you have the dot operator for element-wise operations in play. That's the problem w/ global, the context in the code elsewhere is difficult to debug.

Use the debugger...

Code stylistic note:

dx=[dx(1);dx(2)]; 

doesn't do anything useful...you've already created the column vector orientation for dx in the preallocation step.

[Jonah Foley]

noted, about the stylistic note, honestly forgot I had that preallocation in there and aded the bit you highlighted just recently! The only other code I have for this project is, firstly the pwm function:

function v=pwm(t,T,d)
v=zeros(1,numel(t)); 
for i=1:numel(t)
k = floor(t(i)./T);
if ((k*T<=t(i)) & (t(i)<(k+d)*T))
    v(i)=1;
end
if (((k+d)*T<=t(i)) & (t(i)<(k+1)*T)) %#ok<*AND2>
    v(i)=0;
end
end
end 

this one runs fine, I have simulated and seen expected results, and then:

global v_in R C L T 
v_in = 1.5;
R = 680;
C = 200e-6;
L = 180e-6;
T = 1e-3; 

t = linspace(0,0.2,2e6);
pwm(t,T,0.1); 

options=odeset('MaxStep', 1e-5);
[t,x]= ode23(@boost,(0:1e-7:0.2), [0 0], options);

To be honest, im out of my depth here - since the diff eq. depends on the pwm function im unsure about simullating them both simulataneously, so this may be far from correct. The only reason I used globals is that some of the code was provided for me, and those variables had been already defined as globals - so I assume its necessary for whatever reason.

Hope this sheds some more light on the original problem though

thank you

Answer 1

[Star Strider]

First, pass the extra parameters as arguments. It is best never to use global variables:

function dx=boost(t, x, v_in, R, C, L, T, d)
dx=zeros(2,1);
if pwm(t,T,d)==1
    dx(1) = v_in./L ;
    dx(2) = -x(2)./R*C ;
end
if pwm(t,T,d)==0 & x(1)>=0
    dx(1) = (v_in - x(2))./L ;
    dx(2) = x(1)./C - x(2)./R*C ;
end
if pwm(t,T,d)==0 & x(1)<0           %#ok<*AND2>
    dx(1) = 0 ;
    dx(2) = -x(2)./R*C ;
end
dx=[dx(1);dx(2)]; 
end

Then in your ode23 call, it becomes:

[T,X] = ode23(@(t,x) boost(t, x, v_in, R, C, L, T, d), tspan, x0);

Second, at least one of the values in:

dx(1) = (v_in - x(2))./L;

could be empty, or a vector, either of which would throw that error.

[Jonah Foley]

You may see my other comment in this post to see the full code, but V_in and L are defined with numerical values, so definitely not vectors, and x(2), as far as im aware, is just the component of a column vector.

x is meant to be a 2 element column vector with components x(1) = i_L and x(2) = v_out, of my circuit. MATLAB doesnt need to know that though.

It is also peculiar since all of the variables x(2), v_in and L are used previously to line 9 and they dont cause any of the same errors

[Star Strider]

When I try to run your code (with my verson of ‘boost’):

Unrecognized function or variable 'd'.

What’s ‘d’?

[Jonah Foley]

d is actually an input to the pwm(t,T,d) function and it is a merely a parameter. If you know about the concept of pwm, then d is the duty cycle of the pwm signal.

When it comes to simulation d is then defined as a value in the range 0 < d < 1

[Star Strider]

If you want to use ‘d’ in your code, it has to exist somewhere in your workspace.

This runs without error:

function v=pwm(t,T,d)
v=zeros(1,numel(t)); 
for i=1:numel(t)
k = floor(t(i)./T);
if ((k*T<=t(i)) & (t(i)<(k+d)*T))
    v(i)=1;
end
if (((k+d)*T<=t(i)) & (t(i)<(k+1)*T)) %#ok<*AND2>
    v(i)=0;
end
end
end 

function dx=boost(t, x, v_in, R, C, L, T, d)
dx=zeros(2,1);
if pwm(t,T,d)==1
    dx(1) = v_in./L ;
    dx(2) = -x(2)./R*C ;
end
if pwm(t,T,d)==0 & x(1)>=0
    dx(1) = (v_in - x(2))./L ;
    dx(2) = x(1)./C - x(2)./R*C ;
end
if pwm(t,T,d)==0 & x(1)<0           %#ok<*AND2>
    dx(1) = 0 ;
    dx(2) = -x(2)./R*C ;
end
dx=[dx(1);dx(2)]; 
end

v_in = 1.5;
R = 680;
C = 200e-6;
L = 180e-6;
T = 1e-3; 
d = rand;
t = linspace(0,0.2,2e6);
pwm(t,T,0.1); 
options=odeset('MaxStep', 1e-5);
tspan = (0:1e-7:0.2);
x0 = [0 0];
[tv,xm] = ode23(@(t,x) boost(t, x, v_in, R, C, L, T, d), tspan, x0);

figure
plot(tv,xm)
grid

Here I defined ‘d’ as a uniformly-distributed random number, since that is by defnition (0<d<1), then pass it as a parameter to ‘boost’ that then passes it as a parameter to ‘pwm’. (I am familiar with duty cycles and to an extent with pulse-width modulation, however ‘d’ entirely escaped both definition and documentation.) If ‘d’ is supposed to change in response to something, that needs to be defined as a function somewhere.

[Jonah Foley]

ok, that's really helpful, thank you very much. d doesnt actually change at all, it really only controls the proportion of how long the switch is on and how long the switch is off in my circuit (which remains constant) - so everything you've done here is exactly how I intended to have it working originally.

Using your corrections it finnaly works, ty.

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

[Jonah Foley]

I have accepted your answer but I have one further question, because im not entirely sure whats being plotted . As far as I understand x is a 2x1 coloumn vector, so when you plot the graph of tv against xm, how does it plot the x vector against time? Ideally id like to plot graphs of each component of the x vector against time.

[Star Strider]

Thank you!

In your code, ‘xm’ is a (2E+6 x 2) matrix. (It is always (N x V) where ‘N’ is the number of evaluations returned and ‘V’ is the row length of the vector of ordinary differential equations, here 2.) The ‘tv’ output is a (2E+6 x 1) vector, or generically (N x 1).

To plot them both in the same axes:

figure
plot(tv,xm)
grid
legend('x_1', 'x_2', 'Location','E')

Otherwise ‘x(1)’ is ‘x(:,1)’, and ‘X(2)’ is ‘x(:,2)’, if you want to plot them in separate axes. The ‘tv’ vector is the same for both.

In this run (noting that ‘d’ is random) the plot would be:

[Jonah Foley]

Ahh, that definitely clears it up, So whats being plotted here is infact what I want, which in this case actually refers to the voltage of the circuit (red) and the current (blue). Thanks again.

[Star Strider]

As always, my pleasure!

Considering that the two variables have completely different units, another option for the plot is:

figure
yyaxis left
plot(tv, xm(:,1))
xlabel('Time')
ylabel('Amperes')

yyaxis right
plot(tv, xm(:,2))
ylabel('Volts')

legend('Current', 'Potential', 'Location','NW')

Change the xlabel and ylabel arguments to reflect the actual uints, and the legend arguments to reflect the correct variable names.

See the documentation on yyaxis for details on other changes you can make to the yyaxis plot.

EDIT — (4 Jan 2020 at 19:12)

Added yyaxis plot image —

[Jonah Foley]

Damn, thats looks good thanks!

[Star Strider]

As always, my pleasure!