I don't want any imaginary number , and I don't want any negative number in my matrix, but none of MathWork can even answer it.

I tried to generate the calculation of a (2x10001x5) matrix, with a for loops,

but it continuously gave me the negative and the complex number that I am not asking for.

• if i used the isreal(T) ==0 equation, the entire number becomes zero and the calculation becomes inaccurate, the complex becomes zero and I missed half of the calculations
• I cannot turn the negative number into zero with T(T<0) = 0; I need to eliminate every single negative number whenever the loop goes from w=1:10000. I have no problem if I want to prevent one row or one column from becoming negative or imaginary number, but I can't do the same thing if i want to automatically turn negative or complex number (1,w,:) into zero.

Please let me know how I can incorporate the function that will detect every negative number and turn it into zero, within the loop.

The following is the code

    
%---------setup----------%
x = 4; %numbers of reservoir
S0 = 1.438; % initial storage (mm)
dt = 1; %Timestep (day)
R0 = 2; %precipitation = 2mm/day, it was selected because it was the aveage rainfall in Minnesota
A = 1.202/(1.438*24); %Linear reservoir coefficent (1/day) 
Q0 = S0*A;  %initial outflow (mm/day)
y=1;

A1= rand(x+1,1);
A1 = A1*A;
S1= rand(x+1,1);
S1 = S1*S0;
Q1=rand(x+1,1);
Q1=Q1*Q0;


t=20000
nu01=zeros(2,t/dt+1,x+1);
nu01(1,1,:)=S1; %initial storage when t=1
nu01(2,1,:)=Q1; %initial outflow when t=1
y=1;
Vmax = 7.9022
up=Vmax; % Upper limit
lo = 0; % Lower Limit

%---------setup ends----------%

%The loop I need your help to review
for w = 1:t/dt %days of the experiment; time increment
    nu01(1,w+1,1) =nu01(1,w,1)+Vmax*dt*y-nu01(2,w,1)*dt;                     %storage of Reservoir 1    
    nu01(2,w+1,1) =A1(1)*nu01(1,w+1,1).^N1(1);                                   %outflow of Reservoir 1
    for z = 2:x+1
        nu01(1,w+1,z) =nu01(1,w,z)+Vmax*dt*y-nu01(2,w,z)*dt+nu01(2,w,z-1)*dt;      %storage of Reservoir 2
        nu01(2,w+1,z) =A1(z)*nu01(1,w+1,z).^N1(z);   %outflow of Reservoir 2
    end
    if nu01(2,w,z)>(up-(1e-5)) %after reaching maximum storage; precipitation was turned off
        y=0;
    end
    if sum(nu01(1,w,:)<0)>0 %after reaching minimum storage
        nu02=nu01(1,w,:);
        nu02(nu02<0)=0;
        nu01(1,w,:)=nu02;
    end
    if nu01(2,w,end) < (lo+(1e-5))
        break
    end
end

[John D'Errico]

It is not that nobody can answer your question. In fact, the answer seems pretty simple.

However, we cannot help you at this point. Why not? You don't tell us what is in N1. No place in the code you show is N1 ever defined. So we cannot test your code. However, it is the value of N1 that is almost certainly the crux of the matter.

So it you want help, then you need to make it possible to be helped. What is N1?

Answer 1

[Steven Lord]

I have not read through all your code, but I think there are a number of functions or pieces of functionality that may be of interest or use to you.

• the real and/or max functions
• logical indexing to replace negative or complex values with other values.
• the realpow and possibly realsqrt functions

rng default % So I know exactly which random numbers I'm going to get
sampleData = randi([-10 10], 1, 5)

sampleData = 1×5

7 9 -8 9 3

Let's perform the square root computation then take just the real part of the result.

s1 = sqrt(sampleData)

s1 = 
   2.6458 + 0.0000i   3.0000 + 0.0000i   0.0000 + 2.8284i   3.0000 + 0.0000i   1.7321 + 0.0000i

s1 = real(s1)

s1 = 1×5

2.6458 3.0000 0 3.0000 1.7321

Let's replace the complex values in the result of the square root computation with a placeholder. I'll use NaN.

s2 = sqrt(sampleData)

s2 = 
   2.6458 + 0.0000i   3.0000 + 0.0000i   0.0000 + 2.8284i   3.0000 + 0.0000i   1.7321 + 0.0000i

s2(imag(s2) ~= 0) = NaN

s2 = 1×5

2.6458 3.0000 NaN 3.0000 1.7321

Let's replace the negative numbers in the data with 0 (in two different ways) and then compute using the modified data. I'll leave the original sampleData unchanged by making copies of it.

sampleData2 = sampleData;
sampleData2 = max(sampleData2, 0)

sampleData2 = 1×5

7 9 0 9 3

s3 = sqrt(sampleData2)

s3 = 1×5

2.6458 3.0000 0 3.0000 1.7321

sampleData3 = sampleData;
sampleData3(sampleData3 < 0) = 0

sampleData3 = 1×5

7 9 0 9 3

s4 = sqrt(sampleData3)

s4 = 1×5

2.6458 3.0000 0 3.0000 1.7321

Answer 2

[Bruno Luong]

Just add the max statement to clamp the variable nu01 whenever you set nu01 (2 places):

nu01(1,w+1,z)  = ...;
nu01(1,w+1,z)  = max(nu01(1,w+1,z), 0); % clamp

The test later is redundant (and too late to fix the propagation of negative value) so your for-loop becomes

%The loop I need your help to review
for w = 1:t/dt %days of the experiment; time increment
    nu01(1,w+1,1) = nu01(1,w,1)+Vmax*dt*y-nu01(2,w,1)*dt;                     %storage of Reservoir 1    
    nu01(1,w+1,1) = max(nu01(1,w+1,1), 0); % clamp
    nu01(2,w+1,1) = A1(1)*nu01(1,w+1,1).^N1(1);                                   %outflow of Reservoir 1
    for z = 2:x+1
        nu01(1,w+1,z) = nu01(1,w,z)+Vmax*dt*y-nu01(2,w,z)*dt+nu01(2,w,z-1)*dt;      %storage of Reservoir 2
        nu01(1,w+1,z) = max(nu01(1,w+1,z), 0); % clamp
        nu01(2,w+1,z) = A1(z)*nu01(1,w+1,z).^N1(z);   %outflow of Reservoir 2
    end
    if nu01(2,w,z)>(up-(1e-5)) %after reaching maximum storage; precipitation was turned off
        y=0;
    end
    if nu01(2,w,end) < (lo+(1e-5))
        break
    end
end

[dpb]

Great job of not letting the forest get in the way of the trees in seeing an almost trivial way to solve the problem, Bruno!

Answer 3

[dpb]

    
%---------setup----------%
x = 4; %numbers of reservoir
S0 = 1.438; % initial storage (mm)
dt = 1; %Timestep (day)
R0 = 2; %precipitation = 2mm/day, it was selected because it was the aveage rainfall in Minnesota
A = 1.202/(1.438*24); %Linear reservoir coefficent (1/day) 
Q0 = S0*A;  %initial outflow (mm/day)
y=1;

A1= rand(x+1,1);
A1 = A1*A;
S1= rand(x+1,1);
S1 = S1*S0;
Q1=rand(x+1,1);
Q1=Q1*Q0;


t=20000;
nu01=zeros(2,t/dt+1,x+1);
nu01(1,1,:)=S1; %initial storage when t=1
nu01(2,1,:)=Q1; %initial outflow when t=1
y=1;
Vmax = 7.9022

Vmax = 7.9022

up=Vmax; % Upper limit
lo = 0; % Lower Limit

%---------setup ends----------%

%The loop I need your help to review
for w = 1:t/dt %days of the experiment; time increment
    nu01(1,w+1,1) =nu01(1,w,1)+Vmax*dt*y-nu01(2,w,1)*dt;                     %storage of Reservoir 1    
    nu01(2,w+1,1) =A1(1)*nu01(1,w+1,1).^N1(1);                                   %outflow of Reservoir 1
    for z = 2:x+1
        nu01(1,w+1,z) =nu01(1,w,z)+Vmax*dt*y-nu01(2,w,z)*dt+nu01(2,w,z-1)*dt;      %storage of Reservoir 2
        nu01(2,w+1,z) =A1(z)*nu01(1,w+1,z).^N1(z);   %outflow of Reservoir 2
    end
    if nu01(2,w,z)>(up-(1e-5)) %after reaching maximum storage; precipitation was turned off
        y=0;
    end
    if sum(nu01(1,w,:)<0)>0 %after reaching minimum storage
        nu02=nu01(1,w,:);
        nu02(nu02<0)=0;
        nu01(1,w,:)=nu02;
    end
    if nu01(2,w,end) < (lo+(1e-5))
        break
    end
end

Unrecognized function or variable 'N1'.

We can't run your code to see what you're talking about explicitly; the exponent term in

nu01(1,w+1,z) =nu01(1,w,z)+Vmax*dt*y-nu01(2,w,z)*dt+nu01(2,w,z-1)*dt;   %storage of Reservoir 2
nu01(2,w+1,z) =A1(z)*nu01(1,w+1,z).^N1(z);                              %outflow of Reservoir 2

isn't defined.

BUT, I'm guessing the above is where your problems start; if/when nu01 goes negative and N1 is a fractional power, then nu01^N1 will be complex.

The far better fix is to check for having emptied the reservoir first such that you have a fixup code that prevents the negative to begin with; whether it's adequate for model purposes to just then set the amount to zero or whether you need to add sufficient logic to determine what the actual time between iteration time steps is that matches the zero and keep that additional time point is up to you and how much complexity you want to add.

NOTA BENE:

>> nu01=rand(4,1)
nu01 =
    0.8147
    0.9058
    0.1270
    0.9134
>> nu01(1)=sqrt(-3)
nu01 =
   0.0000 + 1.7321i
   0.9058 + 0.0000i
   0.1270 + 0.0000i
   0.9134 + 0.0000i
>> 

That, if you don't catch the change in sign when it happens and take action then to avoid applying the exponent to a negative term, that when you do and generate that one complex value, MATLAB has to convert the whole array to complex because MATLAB numeric arrays have to be one data type and one datatype only.

But, the answer is to NOT to let it generate the complex numbers in the first place, not to try to patch them afterwards.