I am currently trying to run a numerical model of Daisyworld (which should in theory be quite a simple model). I am using the values from the original Lovelock paper. I am trying to solve it using the forward euler method upon the following equations:
(dαW)/dt=αw (αg β(Tw )-γ)
(dαb)/dt=αb (αg β(Tb )-γ)
I am unsure if I have set out my equations correctly and I am receiving the error :
Here is my code:
clear% clear the workspace
%set initial conditions and constants
q = 2.06e9; %heat flux coefficient k-4
white_albedo = 0.75; %set white albedo
black_albedo = 0.25; %set black albedo
Albedo_bare = 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
P = 1; %proportion of fertile ground set to unity
%preallocate arrays for changing variables
dfw_dt= zeros(1,length(t)) %preallocate array
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
T4e = zeros(length(t)) %preallocate array
Albedo_planetary = zeros(1,length(t)) %preallocate array
%set initial valies of changing variables
dfw_dt(1) = 0
white_area(1) = 0.01
black_area(1) = 0.01
bare_area(1) = 0.98
T4w(1) = 0
BrW(1) = 0
T4e(1) = 0
Albedo_planetary(1) = 0.98*0.5
for i = 2:length(t)
white_area(i+1) = white_area(i) + dfw_dt(i) * tstep
dfw_dt(i) = white_area(i) * bare_area(i) * BrW(i) * T4w(i) - DRate
BrW(i) = 1 - k * (T4w(i) - TempOpt)^2
if (T4w < Topt) BrW = 0
end
T4w(i) = (q * (Albedo_planetary(i) - white_albedo) + T4e(i)^4)
T4e(i) = Solar_lumin * L * (1-Albedo_planetary(i))
Albedo_planetary(i) = (white_area(i)*white_albedo) + (black_area(i) * black_albedo) + (bare_area(i) * Albedo_bare)
bare_area(i) = P - black_area(i) - white_area(i)
end
Any help would be much appreciated!
Thanks
