MATLAB Answers

Complex numbers from for loop (daisyworld model)

2 views (last 30 days)
Jake Lloyd Newman
Jake Lloyd Newman on 16 Feb 2021
Commented: Jon on 16 Feb 2021
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)

Accepted Answer

Jon
Jon on 16 Feb 2021
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.
  2 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by