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

5 次查看(过去 30 天)
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
  8 个评论
Ameer Hamza
Ameer Hamza 2020-3-20
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
Christian Gärtner 2020-3-23
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.

请先登录,再进行评论。

回答(1 个)

Ameer Hamza
Ameer Hamza 2020-3-24
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
  1 个评论
Christian Gärtner
Christian Gärtner 2020-3-26
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!

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Calculus 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by