ode15s integrating over arbitrary function
Show older comments
Hi there,
Im trying to solve this very simple differential equation.
y(t)'=r(T(t))
where r is a rate curve (development as a function of temperature) and T is an arbitrary Temperature profile. I am trying to solve this but i get this error.
%-------------------------------------------------------%
Error using griddedInterpolant
The grid vectors are not strictly monotonic increasing.
Error in interp1 (line 191)
F = griddedInterpolant(X,V,method);
Error in test_ode_mass>Mass (line 97)
M=interp1(x(1,:),x(2,:),t);
Error in test_ode_mass>oderatecurveVar (line 90)
dy(1)=Mass(t,x);
Error in odearguments (line 88)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode15s (line 149)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in test_ode_mass (line 59)
[a cumulativeDevelopment(i,:)] = ode15s(@oderatecurveVar,t,ynot(i), odeset('RelTol',1e-10),
[rT(i,:); t]);
%-------------------------------------------------------%
I have attached to code to the bottom. The main goal is to find a way to solve the differential equation for an arbitrary Temperature profile only specified by a set of data points (not a function). ie a vector of data points relating temperature to time.
%-------------------------------------------------------%
%test new ode solving method
function test_ode_mass
clear all
format long
%============================================================================ %
% time in days
td=.05;
tend=365;
t=0:td:tend;
% temperature variations
T0=10; %average yearly temperature
T1=13; %average fluxuation in temperature per year
Tf=0; %average fluxuation in temperature per day
tPerturb=0; % perturb the temprature by this amount;
% % temperature variations
% T0=5; %average yearly temperature
% T1=15; %average fluxuation in temperature per year
% Tf=8; %average fluxuation in temperature per day
% tPerturb=0; % perturb the temprature by this amount;
% temperature spike variables
pertTempAmp=10;
pertTempExpTime=50;
pertTempDev=0.7;
pertTemperatur=pertTempAmp*exp(-(((t-pertTempExpTime)/pertTempDev).^2)/2);
% perturb the temprature by this amount;
tPhaseShift=pi; % this parameter has no effect on the Gfunction2
tPhaseShift=round(tPhaseShift*10^7)/10^7;
% Temperature spike
T=(T0-T1*cos(2*pi*t/365-tPhaseShift)+(Tf)*cos(pi*t*2)+tPerturb+pertTemperatur);
%TEMPERATURE PLOT
figure
% subplot(3,1,1)
plot(t,T)
title('Temperature as a function of time')
xlabel('Time (days)')
ylabel('Temperature (K)')
axis([min(t) 365 min(T) max(T)])
xr(1,:)={67,8,1.6,1.7,20, T};
xr(2,:)={100,7,1.4,1.8,20.0,T};
xr(3,:)={66,5.5,2,0.4,2.4,T};
xr(4,:)={3,3.5,2.5,0.5,3.5,T};
ynot=zeros(1,length(xr(:,1)));
% options=odeset('RelTol',1e-6,'Stats','on')
for i=1:length(xr(:,1))
rT(i,:)=findRateCurveVar(xr(i,:),t);
[a cumulativeDevelopment(i,:)] = ode15s(@oderatecurveVar,t,ynot(i), odeset('RelTol',1e-10), [rT(i,:); t]);
end
end
%--------------------------------------------%
function rT=findRateCurveVar(x,t)
%%Find Theoretical Developmental rate curves YURK POWELL 2009
%
% Assume that the funcitonal form is
%
% a; maximum developmental time
% b; changes the width of the deveoplental times
% c; slop at lower tempeature threshold
% d; slope at higher threshold temp
% theta; Optimal gorwing temperature (@ minimum develpmental time)
a=x{1};
b=x{2};
c=x{3};
d=x{4};
theta = x{5};
T=x{6};
rT=(1./(a+exp(b-theta*(d+c)+d.*(T)))).*heaviside(T-theta)+(1./(a+exp(b-c.*(T)))).*heaviside(theta-T);
end
%--------------------------------------------%
function dy = oderatecurveVar(t,y,x)
dy=zeros(1,1);
dy(1)=Mass(t,x);
end
%--------------------------------------------%
function M=Mass(t,x)
M=interp1(x(1,:),x(2,:),t);
end
Answers (0)
Categories
Find more on Programming 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!