How to loop through a function twice?

I have a function that now works for individual values. I am trying (and so far, failing) to run the function for several values for each of two parameters: K and Smax.

Any help would be very appreciated. This is my first time using MatLab.

Here is the functional, version I'm using for individual trials:

No Spill Scenario Storage where k=0.3; Smax=100

SynthK=0.3;
SynthSmax=100;
Storage = [Precipmmday(1) 0];
Baseflow = [Precipmmday(1) 0];
for k1 = 2:length(Precipmmday)
    Baseflow (k1,:) = [(Storage(k1-1,2)+Precipmmday(k1-1))*SynthK];
    Storage(k1,:) = [(Storage(k1-1,2)+Precipmmday(k1-1))*(1-SynthK)];    
end
SynthStorage = Storage;
SynthStorage (:,1) =[];

SynthSpill=SynthStorage - SynthSmax;
   SynthSpill(SynthSpill<0)=0;

SynthStorage(SynthStorage>SynthSmax)=SynthSmax;

SynthBaseflow = Baseflow;
SynthBaseflow (:,1) = [];

SynthOutflow = SynthBaseflow + SynthSpill;

And here is what I'm working on to run through many values for the two parameters:

Run Model

  Storage = [Precipmmday(1) 0];
  Baseflow = [Precipmmday(1) 0];
  for Smax = linspace(10,500,10);
      Smax = repmat(Smax,10,1);
      for K = linspace(0.01,1,10);
          K = repmat(K,10,1);
          for k1 = 2:length(Precipmmday)
                  Baseflow(k1,:)  = [(Storage(k1-1,2)+Precipmmday(k1-1))*K];
                  Storage(k1,:) = [(Storage(k1-1,2)+Precipmmday(k1-1))*(1-K)];
          end
              ModelStorage = Storage;

                ModelSpill=ModelStorage - Smax;
                ModelSpill(ModelSpill<0)=0;

                ModelStorage(ModelStorage>Smax)=Smax;

                ModelBaseflow = Baseflow;

                ModelOutflow = ModelBaseflow + ModelSpill;
        end
    end

I can see where some of the errors are but don't know how to fix them. I think the way I'm defining Storage and Baseflow is an issue. Also, I keep getting the "Subscripted assignment dimension mismatch." or "Attempted to access Storage(2,2); index out of bounds becausesize(Storage)=[10,1]." error messages.

Thank you for any suggestions you have.

ETA: Data file is attached

Answer 1

[Star Strider]

You needed to create separate loops for ‘Smax’ and ‘K’. This runs. I will let you decide if it gives the results you want:

Storage = [Precipmmday(1) 0];
Baseflow = [Precipmmday(1) 0];
Smax = linspace(10,500,10);
K = linspace(0.01,1,10);
for k3 = 1:length(Smax)
    for k2 = 1:length(K)
        for k1 = 2:length(Precipmmday)
            Baseflow(k1,k2,k3)  = [(Storage(k1-1,2)+Precipmmday(k1-1))*K(k2)];
            Storage(k1,k2,k3) = [(Storage(k1-1,2)+Precipmmday(k1-1))*(1-K(k2))];
        end
        ModelStorage = Storage;

          ModelSpill=ModelStorage - Smax(k3);
          ModelSpill(ModelSpill<0)=0;

          ModelStorage(ModelStorage>Smax(k3))=Smax(k3);

          ModelBaseflow = Baseflow;

          ModelOutflow = ModelBaseflow + ModelSpill;
      end
  end

I still had the data file and code from yesterday, so I used it to compute ‘Storage’. I tested these revisions on your data.

[Vert]

Thank you!!

[Star Strider]

My pleasure!

I also learned something about surface water hydrology.

[Vert]

Something everyone should know about! :)

I have another question, maybe I should post it separately but it is for the same data.

I am trying to show the level of error as both K and Smax change but right now, only the variation in Smax is shown in the figure I'm generating. Can you see where in the code I am creating this problem?

No Spill Scenario Storage where k=0.3; Smax=100

SynthK=0.3;
SynthSmax=100;
Storage = [Precipmmday(1) 0];
Baseflow = [Precipmmday(1) 0];
for k1 = 2:length(Precipmmday)
    Baseflow (k1,:) = [(Storage(k1-1,2)+Precipmmday(k1-1))*SynthK];
    Storage(k1,:) = [(Storage(k1-1,2)+Precipmmday(k1-1))*(1-SynthK)];    
end
SynthStorage = Storage;
SynthStorage (:,1) =[];

SynthSpill=SynthStorage - SynthSmax;
   SynthSpill(SynthSpill<0)=0;

SynthStorage(SynthStorage>SynthSmax)=SynthSmax;

SynthBaseflow = Baseflow;
SynthBaseflow (:,1) = [];

SynthOutflow = SynthBaseflow + SynthSpill;

Run Model (10x10)

Baseflow = zeros(10,10);
Storage = zeros(10,10);
Storage = [Precipmmday(1) 0];
Baseflow = [Precipmmday(1) 0];
Smax = linspace(10,500,10);
K = linspace(0.01,1,10);
for k3 = 1:length(Smax)
    for k2 = 1:length(K)
        for k1 = 2:length(Precipmmday)
            Baseflow(k1,k2,k3)  = [(Storage(k1-1,2)+Precipmmday(k1-1))*K(k2)];
            Storage(k1,k2,k3) = [(Storage(k1-1,2)+Precipmmday(k1-1))*(1-K(k2))];
        end
        ModelStorage = Storage;
          ModelSpill=ModelStorage - Smax(k3);
          ModelSpill(ModelSpill<0)=0;
          ModelStorage(ModelStorage>Smax(k3))=Smax(k3);
          ModelBaseflow = Baseflow;
          ModelOutflow = ModelBaseflow + ModelSpill;
      end
end

Calculate Difference

SynthOutflow = repmat(SynthOutflow,[1 10 10]);
Difference = ModelOutflow - SynthOutflow;

Root Mean Square Error

rmse = sqrt((sum((Difference).^2))/1087);
rmse1 = reshape(rmse,10,10);

RMSE Figure

figure;
mesh(K, Smax, rmse1);
xlabel('k')
ylabel('Smax')
zlabel('err')

[Star Strider]

I agree!

You need to change this line to limit it to the first 1087 elements, at least because you’re dividing it by that number:

rmse = sqrt((sum((Difference(1:1087,:,:)).^2))/1087);

If you use the entire ‘Difference’ matrix, you get a NaN result because some of the elements of the matrix >1087 are NaN. Make that change and everything works.

[Vert]

Thanks!

I'm still getting the same figure, it shows the difference as SMax changes but nothing is changing along the K axis.

[Star Strider]

My pleasure!

I don’t sufficiently understand what you are doing to suggest how you might change your code. Is it possible that ‘K’ has no effect, or that it is not noticeable with respect to ‘Smax’?

[Vert]

No, I know it will have an effect. I'll keep trying.

[Star Strider]

See my answer to your subsequent Question to follow up on this.