# for loop with stored variables

1 view (last 30 days)
Kirsty James on 14 Mar 2013
Hi i have written a for loop for several variables that i have defined in other coding. I cannot seem to get matlab to recognise these vectors that i have already defined, I am relatively new to matlab and i am unsure if i have set up the for loop correctly. Any help would be great
this defines daisyworld 1
function dydt =daisyworld1(t,y,Lt,L)
Lint=interp1(Lt,L,t); %interpolate the data set (Lt,L) at time t
dydt = [0;0];
A=((1-y(1)-y(2))*0.5)+(y(1)*0.25)+(y(2)*0.75); %albedo
S=917; %constant solar energy
Z=5.67*10^(-8); %Stefan-Boltzmann constant;
Te=((((S*Lint)/Z)*(1-A)).^(1/4))-273; %plantary temperature
B= 1-0.003265*((22.5-Te).^2); %beta value,
g=0.3; %death term
dydt(1)=y(1)*(( (1-y(1)-y(2)) *B)-g); %black daisy formula
dydt(2)=y(2)*(( (1-y(1)-y(2)) *B)-g); %white daisy formula
I then solve it using ode45
clear; % Remove stored variables daisyworldode45
Lt=linspace(0,1000,25); %generate t for L
L=Lt/500+ 0.4; %luminosity function want to keep within a range of 0 and 2 so change depending on tspan
tspan=[0 1000]; %solve for values of t
IC = [0.2 0.3]; %initial conditions of daisy percentage, black then white
[T, Y]=ode45(@(t,y)daisyworld1(t,y,Lt,L),tspan,IC); % solves equation, need capital T and Y
I am trying to extract these values of T and Y, and use the vector T to find values for L and then plot these against Te
N=7269;
Te= zeros(1,N);
for T=0:N
for Y=0:N
Te(T,Y)=((917*(T/500+0.4)/(5.67*10^(-8))*(1-(((1-y(1)-y(2))*0.5)+(y(1)*0.25)+(y(2)*0.75))).^(1/4)))-273;
end
end
plot (T,Y);
7269 is the number of y values that ode45 uses i keep getting the error code
Undefined function 'y' for input arguments of type 'double'.
Error in daisyworldode45 (line 23)
Te(T,Y)=((917*(T/500+0.4)/(5.67*10^(-8))*(1-(((1-y(1)-y(2))*0.5)+(y(1)*0.25)+(y(2)*0.75))).^(1/4)))-273;
Jan on 14 Mar 2013
Remark: 5.67*10^(-8) requires one multiplication and an expensive power function. 5.67e-8 is parsed once during the reading of the M-file and in consequence much more efficient - and nicer.

ChristianW on 14 Mar 2013
There are several ways. I would make for the Plantary Temperature (Te) a new function, and vectorize Albedo (A). Like this:
function test
clear % Remove stored variables daisyworldode45
tspan = [0 1000]; %solve for values of t
L = @(t) t/range(tspan)*2 + 0.4; %luminosity function want to keep within a range of 0 and 2 so change depending on tspan
IC = [0.2 0.3]; %initial conditions of daisy percentage, black then white
[T,Y] = ode45(@(t,y)daisyworld1(t,y,L),tspan,IC); % solves equation, need capital T and Y
LT = L(T);
Te = Temp(LT,Y);
subplot(211), plot(T,Y), xlabel('T, Time'), legend('Y_1','Y_2')
subplot(212), plot(LT,Te)
xlabel('L(T), Luminosity'), ylabel('Te, Plantary Temperature')
function dydt = daisyworld1(t,y,L)
dydt = [0;0];
Lt = L(t); % calculate L at time t
Te = Temp(Lt,y);
B = 1-0.003265*((22.5-Te).^2); % beta value
g = 0.3; % death term
dydt(1) = y(1)*(( (1-y(1)-y(2)) *B)-g); % black daisy formula
dydt(2) = y(2)*(( (1-y(1)-y(2)) *B)-g); % white daisy formula
function Te = Temp(Lt,y)
if isvector(y), y = y(:)'; end % make row vector
A = ((1-y(:,1)-y(:,2))*0.5)+(y(:,1)*0.25)+(y(:,2)*0.75); % albedo
S = 917; % constant solar energy
Z = 5.67*10^(-8); % Stefan-Boltzmann constant;
Te = ((((S*Lt)/Z).*(1-A)).^(1/4))-273; % plantary temperature
Kirsty James on 14 Mar 2013
Thank you so much this has really been a great help, thank you for taking the time to help me.