ode45 error Unable to perform assignment because the indices on the left side are not compatible with the size of the right side

Hi Guys,

I am trying to solve a set of 3 differential equations using ode45. But i can't get it to work because I am getting the following error:

%Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

%Error in freylaw_stat (line 29)
%dy(1,1)=-Ct+(R+r)*y(3);                 %Mr the constant r is introduced here according to Looft et al.(2018) and is an intended variation of the
%original equation

%Error in Gesamtskript_stat>@(t,y)freylaw_stat(t,y,t_min,Fext_M_P,MVC_Verlauf_P) (line 119)
%[t,y]=ode45(@(t,y) freylaw_stat(t,y,t_min,Fext_M_P,MVC_Verlauf_P),timerange,IC);        %PROZENT MUs

%Error in odearguments (line 90)
%f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

%Error in ode45 (line 115)
  %odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

%Error in Gesamtskript_stat (line 119)
%[t,y]=ode45(@(t,y) freylaw_stat(t,y,t_min,Fext_M_P,MVC_Verlauf_P),timerange,IC);        %PROZENT MUs

following is the code I am using for ode45:

timerange=[T0_min tEnd_min];
 IC=[100; 0; 0]; 
 MVC_Verlauf_P=100;
[t,y]=ode45(@(t,y) freylaw_stat(t,y,t_min,Fext_M_P,MVC_Verlauf_P),timerange,IC);

with Fext_M_P being of the type 741979x1 double and using the following function:

function dy=freylaw_stat(t,y,t_min,Fext_M_P,MVC_Verlauf)
L=10;                                   %Modellspezifischer Faktor
r=15;                                   %rest-recovery-multiplier nach 2018#Looft et al.
F=0.00912;                              %Ermüdungsfaktor nach Frey Law (2012)
R=0.00094;                              %Erholungsfaktor nach Frey Law (2012)

Fext_I=Fext_M_P;

MVC=MVC_Verlauf-y(3);
TL0=(Fext_I/MVC)*100;  
    
if (y(2)<TL0)&(y(1)>(TL0-y(2)))
    Ct=L*(TL0-y(2));                    %if Ma<TL and Mr>(TL-Ma)
elseif (y(2)<TL0)&(y(1)<(TL0-y(2)))
    Ct=L*y(1);                          %if Ma<TL and Mr<(TL-Ma)
else
    Ct=L*(TL0-y(2));                    %if Ma>=TL
end

dy(1,1)=-Ct+(R+r)*y(3);                 %Mr 
dy(2,1)=Ct-F*y(2);                      %Ma
dy(3,1)=F*y(2)-R*y(3);                  %Mf
end

Can someone tell me what I am doing wrong here, and how to change it? Your help is greatly appreciated!

Thank you and best regards,

Christian

[Ameer Hamza]

If you add breakpoint at this line

dy(1,1)=-Ct+(R+r)*y(3);                 %Mr

You will see that the right hand side is a vector of size 741979x1. I guess it is expected to be a scalar.

[Christian Gärtner]

Hi Ameer,

Okay, I understand that.

But actually I want the left hand side and the right hand side to be a vector of size 741979x1. Could you tell me how to do this?

Thanks for taking the time to help me!

[Ameer Hamza]

Can you explain what you are you trying to calculate? I am sure you don't want to make dy of size 741979x1. It will mean that you have a total of 741979 differential equations. If you can describe the actual problem, it will help in suggesting an appropriate solution.

Secondly, even if you have 741979 differential equations, then how are you going to use Ct-F*y(2) and F*y(2)-R*y(3).

[tmarske]

What do you want the function output value dy to be here? You don't explicitly initialise it anywhere, so when it gets to this line:

dy(1,1)=-Ct+(R+r)*y(3);                 %Mr 

Matlab will see the left hand side and assume that dy is a 1*1 numerical array. The RHS is a different size, so Matlab can't assign it to dy(1,1) and throws the error.

What do you want dy to be? A 3 * 741979 matrix? A 3 * 1 cell array of where each cell contains a single 741979x1 vector?

[Christian Gärtner]

Okay, so let me briefly explain: this is in context of human muscle fatigue. Input is the vector 741979x1 vector TLO which is a measure of the force a muscle needs to perform a movement relative to its maximum force output. R,r,F and L are constant fatigue and recovery factors.

As an Output there is

1) Ma which is representing the percentage of force the muscle is actively outputting when it performs a movement. This is calculated by the equation: dy(2,1)=Ct-F*y(2);

2) Mr which is representing the percentage of force of the muscle which is still resting, because it is not needed to perform the movement. This is calculated by the equation: dy(1,1)=-Ct+(R+r)*y(3);

3) Mf s representing the percentage of force of the muscle which is getting fatigued, after it has been active for some time when a movement is performed. This is calculated by the equation: dy(3,1)=F*y(2)-R*y(3);

So, what i want dy to be is probably a 741979 * 3 matrix. Or three single vectors of 741979*1, which represent Ma, Mr and Mf, so I can plot them into a figure.

If this is interesting to you, the equations and further explanations can be found in Xia & Frey Law (2008): A theoretical approach for modeling peripheral muscle fatigue and recovery (https://www.researchgate.net/publication/23255929_A_theoretical_approach_for_modeling_peripheral_muscle_fatigue_and_recovery).

Thank you guys a lot for helping me!

[Christian Gärtner]

Just to give you some more information from the original source, here is a schematic drawing of the function of the muscle fatigue model.

And the original equations I am trying to calculate with dy(1,1), dy(2,1), and dy(3,1):

[Ameer Hamza]

To solve them as a differential equation with ode45, you need dy to be 3*1, not 741979*3. The equation you pasted also show that C(t) is a function of time. So you need to specify, how the time t is related to the elements of vector Ct. For example, if t=0, C(0) is the first element of Ct, i.e., Ct(1), then at the next time step, say t=0.1, which elemnt of Ct is equal to the value of C(0.1).

[Christian Gärtner]

thanks for your answer Ameer!

okay, I did not know that. So I would need to use the following code in the function?

dy(1,:)=-Ct+(R+r)*y(3);                 %Mr the constant r is introduced here according to Looft et al.(2018) and is an intended variation of the original equation
 dy(2,:)=Ct-F*y(2);                      %Ma
 dy(3,:)=F*y(2)-R*y(3);                  %Mf

And for C(t): it has the same length as the time vector so the corresponding value to t=0 is the first element of C. In the next time step the second element of C is the corresponding value, and so on. C(t) is calculated with the following equation:

There LD/LR=10=const, TL is a vector with 741979 elements, where the first cell corresponds to t=0, the second cell corresponds to the next time step and so on. Ma and Mr are the values which are calculated with this differential equation at the corresponding time steps.

I can also send you the whole file with the vectors, if that makes things easier for you.

Answer 1

[Ameer Hamza]

Check the following code. I created a random vector Fext_M_P to test the output of my code. Also the time range is from 0 to 10.

Fext_M_P = rand(741979, 1);
timerange=linspace(0, 10, numel(Fext_M_P));
IC=[100; 0; 0];
MVC_Verlauf_P=100;
[t,y]=ode45(@(t,y) freylaw_stat(t,y,timerange,Fext_M_P,MVC_Verlauf_P),timerange,IC);

function dy=freylaw_stat(t,y,timerange,Fext_M_P,MVC_Verlauf)
L=10;                                   %Modellspezifischer Faktor
r=15;                                   %rest-recovery-multiplier nach 2018#Looft et al.
F=0.00912;                              %Ermüdungsfaktor nach Frey Law (2012)
R=0.00094;                              %Erholungsfaktor nach Frey Law (2012)

index = floor(interp1(timerange([1 end]), [1 numel(Fext_M_P)], t));

Fext_I=Fext_M_P(index);
MVC=MVC_Verlauf-y(3);
TL0=(Fext_I/MVC)*100;

if (y(2)<TL0)&(y(1)>(TL0-y(2)))
    Ct=L*(TL0-y(2));                    %if Ma<TL and Mr>(TL-Ma)
elseif (y(2)<TL0)&(y(1)<(TL0-y(2)))
    Ct=L*y(1);                          %if Ma<TL and Mr<(TL-Ma)
else
    Ct=L*(TL0-y(2));                    %if Ma>=TL
end

dy(1,1)=-Ct+(R+r)*y(3);                 %Mr
dy(2,1)=Ct-F*y(2);                      %Ma
dy(3,1)=F*y(2)-R*y(3);                  %Mf
end

[Christian Gärtner]

Hello there, sorry for the late answer, I am quite busy lately.

Thank you a lot for your answer, it worked great to get a result for the ODE-system!

Still, I was expecting another result, for dy(3,1). It is not behaving as expected, although dy(1,1) and dy(2,1) are perfectly reasonable results.

I have attached a picture of the plotted results where MR is dy(1,1); MA is dy(2,1) and MF is dy(3,1). Also I have drawn a yellow line, as dy(3,1) is expected to behave.

Maybe you have got any idea, what this different behaviour is related to? I will check units and stuff, but I think that they are right.

What I think may be the issue is that MVC is not decreasing over the timerange? Actually it should decrease in the form of:

MVC=MVC_Verlauf-y(3)

But I am unsure how to check MVC, because I do not know how to take a look at it.

Thank you!