MATLAB: Forward euler with multiple changing variables (Daisyworld model)

Question

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

Answer 1

Hi there,

I believe the error you are seeing is has to do with how the variable 'BrW' is being set. While 'BrW' has been initially set as an array, in the following line within the for-loop:

if (T4w < Topt) BrW = 0

the entire of 'BrW' is being updated as a double variable. Then on subsequent iterations of the for-loop, when function call attempts to access 'BrW' as if it were an array of elements, the call fails with the error message you are seeing.

Try replacing the above line with:

if (T4w < TempOpt)  BrW(i) = 0

which tells the loops to only update the specific entry in 'BrW', rather than the entire variable.

Cheers,

Kevin