Including changing variable in ODE

Question

0 个投票

Hi,

I have writted an ODE that worked when I used alpha (Earth's albedo) as a constant, however, now I am trying to vary the albedo with respect to temperature. I wrote a simple equation to vary it in line with what was set out in the question, and it worked...

m=-0.02;
j=linspace(-20,20,100);
alpha_variable = 0.5+m*j; %using bounding albedo and corrected by 0.5 so it goes through horizontal lines
alpha_variable(j<-10) = 0.7; 
alpha_variable(j>10) = 0.3; 
figure(3)
plot(j,alpha_variable)
xlabel('Temperature variations')
ylabel('albedo (alpha)')
title('Albedo response to temperature')

...but now when I try to put in my alpha_variable into the ODE I get an error that the index exceeds the array bounds.

alpha_v=0.7; %I've set alpha_v to 0.7 as this is the value it would be at my lowest Tsvariable
for i=1:9 %9 columns in S_variable
  j=1:length(alpha_v);
  emb_variable=@(t,Ts)((((1-alpha_v(1,j))*S_variable(:,i)/4)-(e_lw*sigma*Ts^4))/c_eff);
   [t,Tsvariable]=ode45(emb_variable,(0:3e7:50*3e7),0);
   Ts_variable=Tsvariable-273;
   m=-0.02;
    alpha_v = 0.5+m*Ts_variable; %using bounding albedo and corrected by 0.5 so it goes through horizontal lines
    alpha_v(Ts_variable<-10) = 0.7; 
    alpha_v(Ts_variable>10) = 0.3;
   figure(2);
   hold on
   plot(t,Ts_variable)
end
xlabel('time (yrs)');
ylabel('temperature (Celcius)');
title('Temperature as a function of time for variable solar constants')
legend('Solar constant = 2734W/m^2','Solar constant = 2460/m^2','Solar constant = 2187W/m^2','Solar constant = 1913W/m^2','Solar constant = 1640W/m^2','Solar constant = 1367W/m^2 (present day)','Solar constant = 1093W/m^2','Solar constant = 820W/m^2','Solar constant = 547W/m^2','Location','southeast')