cannot get my equations to recognise my vector syntax

Hi I am working on some coding and i'm struggling to get the syntax for my vectors correct here are my codes

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.4, 0.3]; %initial conditions of daisy percentage, black then white

[T,Y] = ode45(@(t,y)daisyworld2(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'), ylabel( 'Area'), legend('Black Daisies', 'White Daisies')
subplot(212), plot(LT,Te)
xlabel('L(T), Luminosity'), ylabel('Te, Plantary Temperature')

and another code

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

and the final coding

function dydt = daisyworld2(t,y,L)

dydt = [0;0];

Lt = L(t); % calculate L at time t

Te = Temp(Lt,y);

A = ((1-y(:,1)-y(:,2))*0.5)+(y(:,1)*0.25)+(y(:,2)*0.75); % albedo

Tb= (((2.06*10^9)/(4*((273+22.5).^3)))*(A-0.75))+Te;
Tw= (((2.06*10^9)/(4*((273+22.5).^3)))*(A-0.25))+Te;

Bb = 1-0.003265*((22.5-Tb).^2); % beta value
Bw = 1-0.003265*((22.5-Tw).^2); % beta value

g = 0.3; % death term

dydt(1) = y(:,1)*((  (1-y(:,1)-y(:,2))  *Bb)-g); % black daisy formula
dydt(2) = y(:,2)*((  (1-y(:,1)-y(:,2))  *Bw)-g); % white daisy formula

The error message i keep receiving is

Attempted to access y(:,2); index out of bounds because size(y)=[2,1].

Error in daisyworld2 (line 9)
A = ((1-y(:,1)-y(:,2))*0.5)+(y(:,1)*0.25)+(y(:,2)*0.75); % albedo

Error in test>@(t,y)daisyworld2(t,y,L) (line 10)
[T,Y] = ode45(@(t,y)daisyworld2(t,y,L),tspan,IC); % solves equation, need capital T and Y

Error in odearguments (line 88)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 114)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in test (line 10)
[T,Y] = ode45(@(t,y)daisyworld2(t,y,L),tspan,IC); % solves equation, need capital T and Y

I think its just my vector syntax that i cant get right but any help would be greatly appreciated. Thank you

Answer 1

[Cedric]

In

 A = ((1-y(:,1)-y(:,2))*0.5)+(y(:,1)*0.25)+(y(:,2)*0.75); % albedo

you are trying to access all rows of column 2 of y, but y is a column vector (there is no 2nd column). Without looking too much at what you did, my guess is that you want y(1) and y(2) .. and it is the same below when you define dydt(1) and dydt(2).

[Kirsty James]

I have changed these variables to just y(1) and y(2) but it does not seem to recognise both values of y that i input?

[Cedric]

So you modified also the function Temp, and now you have a result (no more dimension mismatch) which doesn't correspond to what you are expecting?

If so, it might be time to use the debugger and check step by step that what you are computing is correct, or just a few print statements t relevant locations.

[Kirsty James]

i've tried adding the y(:,1) and y(:,2) notation in everytime i use y but i am getting the error message

Error using odearguments (line 91)
@(T,Y)DAISYWORLD2(T,Y,L) must return a column vector.

Error in ode45 (line 114)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in test (line 10)
[T,Y] = ode45(@(t,y)daisyworld2(t,y,L),tspan,IC); % solves equation, need capital T and Y

so i think i am still having trouble with miss matched size

[Kirsty James]

i'm struggling to get the codes to recognise the second entry of my y vector, on the graphs it just plots it as zero no matter what initial condition i enter, any help would be great.

[Cedric]

Did you keep the

 dydt = [0;0];

that you had in DAISYWORLD2? If you dropped it or changed the semi-column for a comma, this might be the reason why it's not working.

The following code, where I performed the changes that I mentioned and only them, should be working actually:

 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.4, 0.3]; %initial conditions of daisy percentage, black then white
    [T,Y] = ode45(@(t,y)daisyworld2(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'), ylabel( 'Area'), legend('Black Daisies', 'White Daisies')
    subplot(212), plot(LT,Te)
    xlabel('L(T), Luminosity'), ylabel('Te, Plantary Temperature')
 end

 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
 end

 function dydt = daisyworld2(t,y,L)
    dydt = [0;0];
    Lt = L(t); % calculate L at time t
    Te = Temp(Lt,y);
    A = ((1-y(1)-y(2))*0.5)+(y(1)*0.25)+(y(2)*0.75); % albedo
    Tb= (((2.06*10^9)/(4*((273+22.5).^3)))*(A-0.75))+Te;
    Tw= (((2.06*10^9)/(4*((273+22.5).^3)))*(A-0.25))+Te;
    Bb = 1-0.003265*((22.5-Tb).^2); % beta value
    Bw = 1-0.003265*((22.5-Tw).^2); % beta value
    g = 0.3; % death term
    dydt(1) = y(1)*((  (1-y(1)-y(2))  *Bb)-g); % black daisy formula
    dydt(2) = y(2)*((  (1-y(1)-y(2))  *Bw)-g); % white daisy formula
 end

[Kirsty James]

I've entered those codes into matlab and i'm still getting a straight line of zero for my white daisy value?

[Cedric]

I get the plots here, and these aren't straight lines. They look quite close to 0 a large part of the time though.

If you place your cursor in function test on the following line

 LT = L(T);

and press F12, it will set a breakpoint there. Then you press F5 to execute your script, which will run until this point. You can then look at Y and for example zoom on the first part of the dynamics by executing

plot(Y(1:100,:))

and you should see that Y starts at the initial condition and decreases. Then you can execute your program step by step pressing F10, or stepping into functions when possible by pressing F11, or going on the executing by pressing F5, or interrupting the debugging session by pressing Shift+F5.