Indexing Error's line 65

3 次查看（过去 30 天）
Zack Bayhan
Zack Bayhan 2014-11-3
评论： Star Strider 2014-11-4
Hello Everyone, I'm going to post the code below but let me explain the issue at hand first. The code will work properly 1 out of 20ish times and will provide a matrix times which should be the moment the concentration of substance b calculated in the with the differential equation. So when the code works properly all 50 passes of the loops achieve the concentration of greater than .6 however the code breaks down when a concentration does not achieve .6 or greater. I attempted to account for this with the while loop at the bottom of the code but I can't quite figure it out. Any tips would be greatly appreciated. Thanks in advance
clear
tic;
% Converting temp to absolute scale and changing to 1/T
T=10:10:100;
T=T+273; T=T'.^-1;
ao=ones(length(T),1);
Z=[ao -T];
% Finding values for k1
k1=[.0859 .131 .243 .396 .646 1.07 1.65 2.50 3.71 5.36];
k1=log(k1);
k1=k1';
% Finding Values for k2
k2=[.0071 .0239 .0718 .205 .553 1.41 3.40 7.81 17.1 36.0];
k2=log(k2);
k2=k2';
%solving for alpha and E/R
x1=Z\k1; % alpha1 and E1
x2=Z\k2; % alpha2 and E2
alp1=x1(1);
alp2=x2(1);
AE1=x1(2);
AE2=x2(2);
% Finding ssr first data set
SSR1=sum((k1-(Z*x1)).^2); %Sum of squared residuals
syx1=sqrt(SSR1/(length(T)-length(x1))); %Standard Error of Estimate
% Finding mean and STD for first data set
f=((inv(Z'*Z)));
var1=f(2,2).*syx1^2;
std1=sqrt(var1);
% Finding ssr 2nd data set
SSR2=sum((k2-(Z*x2)).^2); %Sum of squared residuals
syx2=sqrt(SSR2/(length(T)-length(x2))); %Standard Error of Estimate
% Finding mean and STD for second data set
f=((inv(Z'*Z)));
var2=f(2,2).*syx2^2;
std2=sqrt(var2);
n=51;
for i=1:n
e1(i)=std1*randn(1)+AE1;
e2(i)=std2*randn(1)+AE2;
alpm1(i)=randn(1)+alp1;
alpm2(i)=randn(1)+alp2;
R=8.341;
Tm=1/295;
k1m(i)=alpm1(i)*exp(-e1(i)/R*Tm);
k2m(i)=alpm2(i)*exp(-e2(i)/R*Tm);
% solving the system of differential equations
dcdt=@(t,c)[-k1m(i)*c(1);k1m(i)*c(1)-k2m(i)*c(2)];
[time,c]=ode45(dcdt,[0 10],[1;0]);
conb=c(:,2);
x=0;
while i<50 && x<1
time=find(conb>.6);
times(i)=time(1);
x=1;
if conb(i)<.6 && ~isempty(time)
i=i+1;
x=0;
end
end
end
toc;
  6 个评论
显示 4更早的评论
Image Analyst
Image Analyst 2014-11-4
I repeatedly get the error
Index exceeds matrix dimensions.
Error in test2 (line 65)
times(i)=time(1);
when I try to run it. Like I said, time(1) throws an exception, even if you just type it in on the command line.
Zack Bayhan
Zack Bayhan 2014-11-4
@ Star Strider, That is correct it is a monte carlo type study, the randomaly generated numbers are then used in the differential equation in the for loop. I'm specifically looking for the time that the concentration of substance B increases over .6. When the code works I get a value for all the times in which the vector is greater than .6 which is called time, so the vector times is the first value in which the concentration goes over .6. If this was going to work correctly I should get 50 values for times which I can then do the statistics on. Did this help at all? and thank you for helping with the formatting of the code. I've read the link and hopefully I will be able to format correctly from now on.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2014-11-4
编辑：Star Strider 2014-11-4
I would forget about the threshold in your code and record all the data and parameters for all runs for your simulation in a .mat file. For that to work most efficiently, you would need to define your ode45 ‘tspan’ argument as a vector rather than a range, for example tspan=linspace(0,10,25) to create a 25-element time vector, keeping it the same for all runs. Then after that determine the times ‘conb’ exceeds your threshold.
The advantage to that approach is that your code wouldn’t crash — do perhaps 1000 runs — and you would have all your data at the end to explore at your leisure. That’s how I would do it.
You would have to rewrite parts of your code to do this, but it would likely be worth the effort.
  2 个评论
Zack Bayhan
Zack Bayhan 2014-11-4
alright I think I understand what you are saying however I may have some difficulty getting there. The range you are talking about, is this the bounds of my integration set it the ode45 syntax [0 10]? If I wanted to change this to a linespace instead would I just replace that term in the syntax? and what command would I use to get the data saved in a .mat file?
Star Strider
Star Strider 2014-11-4
It shouldn’t be too difficult. You simply have to save your data from the integration as arrays, and the parameters for each one as elements of vectors, with the subscripts of the ODE solver outputs and the vectors sharing the same subscripts. You can even use cell arrays at every iteration to save all of them together.
If you want to change ‘tspan’ to be a vector of discrete times rather than a range, use linspace:
tspan=linspace(0,10,25);
That gives you 25 points from 0 to 10. If you want a different number, change the 25 to what you want.
With respect to using .mat-files, see the documentation on Save, Load, and Delete Workspace Variables for details.

请先登录，再进行评论。

更多回答（2 个）

Image Analyst
Image Analyst 2014-11-4
It's horrible programming practice to change the iterator variable, "i", inside a for loop. Anyway, what happens is that it changes it for that iteration but then it gets reset to what it should be when that iteration ends. Perhaps the while should use k or some other letter as an iterator - not sure since I don't know the purpose of the while loop due to a glaring lack of comments for it.
  1 个评论
Zack Bayhan
Zack Bayhan 2014-11-4
Sorry about the lack of comments inside the loop. The purpose of the while loop is an attempt to fix a problem I've been running into, before the while loop was in there the the program would function properly assuming that every iteration of the for loop produced a value for conb that was over .6. when the vector for conb didn't produce a value over .6 the code would break down. My attempt to fix this was to carry on the normal operation with the while loop and if the vector failed to produce a value over the .6 figure I wanted the if command to simply add another pass of the loop. I was hoping that once there was 50 values in the times vector the loop would stop.

请先登录，再进行评论。

Zack Bayhan
Zack Bayhan 2014-11-4
Maybe this will be helpful, here is the code prior to the while loop going in.
clear
tic;
% Converting temp to absolute scale and changing to 1/T
T=10:10:100;
T=T+273; T=T'.^-1;
ao=ones(length(T),1);
Z=[ao -T];
% Finding values for k1
k1=[.0859 .131 .243 .396 .646 1.07 1.65 2.50 3.71 5.36];
k1=log(k1);
k1=k1';
% Finding Values for k2
k2=[.0071 .0239 .0718 .205 .553 1.41 3.40 7.81 17.1 36.0];
k2=log(k2);
k2=k2';
%solving for alpha and E/R
x1=Z\k1; % alpha1 and E1
x2=Z\k2; % alpha2 and E2
alp1=x1(1);
alp2=x2(1);
AE1=x1(2);
AE2=x2(2);
% Finding ssr first data set
SSR1=sum((k1-(Z*x1)).^2); %Sum of squared residuals
syx1=sqrt(SSR1/(length(T)-length(x1))); %Standard Error of Estimate
% Finding mean and STD for first data set
f=((inv(Z'*Z)));
var1=f(2,2).*syx1^2;
std1=sqrt(var1);
% Finding ssr 2nd data set
SSR2=sum((k2-(Z*x2)).^2); %Sum of squared residuals
syx2=sqrt(SSR2/(length(T)-length(x2))); %Standard Error of Estimate
% Finding mean and STD for second data set
f=((inv(Z'*Z)));
var2=f(2,2).*syx2^2;
std2=sqrt(var2);
for i=1:50
e1(i)=std1*randn(1)+AE1;
e2(i)=std2*randn(1)+AE2;
alpm1(i)=randn(1)+alp1;
alpm2(i)=randn(1)+alp2;
R=8.341;
Tm=1/295;
k1m(i)=alpm1(i)*exp(-e1(i)/R*Tm);
k2m(i)=alpm2(i)*exp(-e2(i)/R*Tm);
% solving the system of differential equations
dcdt=@(t,c)[-k1m(i)*c(1);k1m(i)*c(1)-k2m(i)*c(2)];
[t,c]=ode45(dcdt,[0 10],[1;0]);
conb=c(:,2);
plot(t,conb)
hold on
time=find(conb>=.6);
times(i)=time(1);
end
toc;
  3 个评论
显示 1更早的评论
Zack Bayhan
Zack Bayhan 2014-11-4
Well that worked out great! Now I have a matrix containing 10000 1x8 matrix's no I just have to figure out how to get matlab to retrieve the desired data. Thanks for all the help.
Star Strider
Star Strider 2014-11-4
My pleasure!
Retrieving the data is relatively straightforward. You can use the load function (always with an output argument so you have some control over the variables you load), but I prefer using the matfile function because it gives me more interactive control over what I load into my workspace. All those have references and links in the link I provided in my Comment, so I won’t repeat them here.
You probably need to do some statistical analyses on your data to determine how best to describe them. With 10000 simulations, going with the normal distribution is probably safe, although it would likely be best to explore some alternatives, such as the lognormal distribution if you know your data are never negative.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by