Clear Filters
Clear Filters

for loop with stored variables

5 views (last 30 days)
Kirsty James
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 Comments
Kirsty James
Kirsty James on 14 Mar 2013
I'm taking the temperature function, Te, from the original daisyworld function. The vectors T and Y are calculated in ode45 and i want to use those values to calculate the Te at each point and plot it against luminosity
Jan
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
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 Comment
Kirsty James
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.

More Answers (0)

Categories

Find more on Mathematics 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!