for loop with stored variables

2 views (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;
3 CommentsShow 2 older commentsHide 2 older comments
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.

Sign in to comment.

Accepted Answer

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
1 CommentShow NoneHide None
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.

Sign in to comment.

Categories

Find more on Logical in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!