Indexing Error's line 65

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;

[Image Analyst]

Read this: http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup, or better yet, just attach your m-file.

How did you get it to go past this line:

times(i)=time(1);

without throwing an exception:

Index exceeds matrix dimensions.
Error in test2 (line 64)
    times(i)=time(1);

time(1) throws an exception even all by itself in the command window.

[Zack Bayhan]

just keep clicking run a bunch of times one out of 20 or so will complete the code. when the elapsed time comes into the command window you will know that the program has finished if it doesn't finish you get the error.

[per isakson]

Delete your comment and learn how to use the {}Code-button. Select text and click the button.

This much code is better to attach with the paper-clip-button.

[Star Strider]

It looks as though you’re doing a Monte Carlo type study, with normally distributed random parameters. You would expect it to not always produce the same results.

What is the significance of the 0.6 threshold? I ask because it might be more interesting from a statistical perspective to store the parameter values and then look at the distribution of conb as functions of them. You can then calculate the probability that conb>0.6.

(I did the reformatting. I’ve seen this code before.)

[Image Analyst]

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]

@ 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.

Answer 1

[Star Strider]

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.

[Zack Bayhan]

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]

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.

Answer 2

[Image Analyst]

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.

[Zack Bayhan]

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.

Answer 3

[Zack Bayhan]

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;

[Star Strider]

I’m not quite sure what you’re doing. To do what I suggested, replace these two lines:

time=find(conb>=.6);
times(i)=time(1);

with:

simulation{i} = {e1(i) e2(i) alpm1(i) alpm2(i) k1m(i) k2m(i) t c};

Then after the end of your (badly-named) ‘i’ loop, save your ‘simulation’ variable (my name for illustration purposes) to a mat-file as:

save('Reaction_Simulations', 'simulation');

Name everything in terms that are meaningful to you. My variable names and file names are simply for illustration.

I really would run several hundred simulations — perhaps as many as 1000 were this my study — and save them. Then you can explore them later without having to re-run your entire code. That will make your life easier, and there is much to be said for that!

See my Comment to my Answer that includes the link to ‘Save, Load, and Delete Workspace Variables’ for details on using the save and load functions.

[Zack Bayhan]

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]

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.