Complex numbers from for loop (daisyworld model)

5 次查看（过去 30 天）
Jake Lloyd Newman
Jake Lloyd Newman 2021-2-16
评论： Jon 2021-2-16
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)

请先登录，再回答此问题。

采纳的回答

Jon
Jon 2021-2-16
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 个评论
Jake Lloyd Newman
Jake Lloyd Newman 2021-2-16
Hi Jon,
Thanks for your help, I was indeed raising negative values to fractional powers!
Cheers,
Jake

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by