Complex numbers from for loop (daisyworld model)

I am currently trying to run a numerical model of Daisyworld. I am using the values from the original Lovelock paper.

When I run my code to calculate growth of the daisies, it returns complex numbers for my variables rather than real numbers. What should I do in order to get real numbers? Any help in solving this would be appreciated.

Here is my code:

clear % clear the workspace
%set initial conditions and constants
q = 2.06*10^9; %heat flux coefficient k-4
white_albedo = 0.75; %set white albedo
black_albedo = 0.25; %set black albedo
bare_albedo = 0.5; %set bare albedo
SB = 5.67e-8; % Stefann Boltzmann constant
Solar_lumin = 917; %solar luminosity (W m-2)
TempOpt = 295.65 ; %this is the optimum temperature for daisy growth in celsius 22.5/ 298
B = 0.25; %constant used in analytical solution
k = 17.5^-2; %width parametetr of daisy temperature growth 290.65K/17.5 celsius
qt= 79.79; % rescaled q
DRate = 0.3; %death rate (gamma)
initial_coverage = 0.01; %set initial coverage of B/W daises as 1% 
black_area = initial_coverage; %initial coverage of black daisies = 0.01
white_area = initial_coverage; %initial coverage of white daisies = 0.01
bare_area = 0.98; %initial coverage of bare ground = 0.98 
L = 1;%set  luminosity 
t_start = 0
t_end = 10
tstep = 1
t=t_start:tstep:t_end; %time series of 0-100
SL = Solar_lumin *L 
P = 1; %proportion of fertile ground set to unity
%preallocate arrays for changing variables
dfw_dt= zeros(1,length(t)); %preallocate array
dfb_dt = zeros(1,length(t));
white_area=zeros(1,length(t)); %preallocate array
black_area=zeros(1,length(t)); %preallocate array
bare_area = zeros(1,length(t)); %preallocate array
T4w = zeros(1,length(t)); %preallocate array
BrW = zeros(1,length(t)); %preallocate array
Tb = zeros(1,length(t));
BrB= zeros(1,length(t));
T4e = zeros(1,length(t)); %preallocate array
Albedo_planetary = zeros(1,length(t)); %preallocate array
%set initial valies of changing variables
%dfw_dt(1) = 0;
%dfb_dt(1) = 0;
white_area(1) = 0.01;
black_area(1) = 0.01;
bare_area(1) = 0.98;
%T4w(1) = 0;
Albedo_planetary(1) = (white_area(1) * white_albedo) + (black_area(1) *black_albedo) + (bare_area(1) * bare_albedo)
T4e(1) = ((SL) / (4*SB ) * (( 1 - Albedo_planetary(1)))) ^0.25
T4w(1) = (q*(Albedo_planetary(1) - white_albedo) + (T4e(1)^4)) ^(0.25)
%Tb(1) = 0
Tb(1) = ((q * (Albedo_planetary(1) - black_albedo) + (T4e(1)^4)))^(0.25)
%BrW(1) = 0;
BrW(1) = (1 - k) * ((T4w(1) - TempOpt))^2
%BrB(1)= 0;
 BrB(1) = ((1- k) * (Tb(1) - TempOpt)) ^2
%T4e(1) = 0;
dfw_dt(1) = white_area(1) * (bare_area(1) * BrW(1) * T4w(1) - DRate)
dfb_dt(1) = black_area(1) * (bare_area(1) * BrB (1) * Tb(1) - DRate)
for i = 2:length(t)
      bare_area(i) = P - black_area(i-1) - white_area(i-1)
      Albedo_planetary(i) = (white_albedo * white_area(i-1)) + (black_albedo * black_area(i-1)) + (bare_albedo * bare_area(i-1))
         T4e(i) = ((SL) / (4*SB ) * (( 1 - Albedo_planetary(i-1)))) ^(0.25)
         T4w(i) = ((q *(Albedo_planetary(i-1) - white_albedo) + ((T4e(i-1)).^4))) ^0.25
          Tb(i) = ((q * (Albedo_planetary(i-1) - black_albedo) + (T4e(i-1).^4)))^0.25
    
    BrW(i) = (1 - k) * ((T4w(i-1) - TempOpt))^2
    if (T4w(i-1) - TempOpt < k^ 0.5) BrW(i) = 0
    end
    BrB(i) = ((1- k) * (Tb(i-1) - TempOpt)) ^2
    if (Tb(i-1) - TempOpt < k^0.5) BrB(i) = 0
    end
    
    dfw_dt(i) = white_area(i-1) * (bare_area(i-1) * BrW(i-1) * T4w(i-1) - DRate)
    dfb_dt(i) = black_area(i-1) * (bare_area(i-1) * BrB (i-1) * Tb(i-1) - DRate)
   
    white_area(i) = white_area(i-1) + dfw_dt(i-1) * tstep
    black_area(i) = black_area(i-1) + dfb_dt(i-1) * tstep
    
end
figure (1)
plot(t,white_area)
hold on
plot(t,black_area)

Answer 1

[Jon]

Without going through your code in detail (there is a lot there) I would suspect that you are raising negative values to fractional powers. For example here where you have an expression raised to the 0.25 power:

(Albedo_planetary(1) - black_albedo) + (T4e(1)^4)))^(0.25).

I would suggest stepping through with the debugger and check the argument values to that and similar expressions to see if they are negative.

If this is just a numerical issue and the imaginary parts are very small compared to the real parts, you could just use real() to keep the real portion of the complex variables that are causing you problems.

[Jake Lloyd Newman]

Hi Jon,

Thanks for your help, I was indeed raising negative values to fractional powers!

Cheers,

Jake

[Jon]

Glad to hear that helped