cannot get my equations to recognise my vector syntax

Question

Kirsty James 2013-3-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/67810-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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Cedric 2013-3-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/67810-cannot-get-my-equations-to-recognise-my-vector-syntax#answer_79219

编辑：Cedric 2013-3-19

在 MATLAB Online 中打开

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).

7 个评论
显示 5更早的评论隐藏 5更早的评论

Cedric 2013-3-19

编辑：Cedric 2013-3-19

在 MATLAB Online 中打开

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 2013-3-19

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

Cedric 2013-3-19

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

cannot get my equations to recognise my vector syntax

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

7 个评论
显示 5更早的评论隐藏 5更早的评论

另请参阅

类别

标签

Community Treasure Hunt

cannot get my equations to recognise my vector syntax

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

7 个评论 显示 5更早的评论隐藏 5更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论