MATLAB: Error using ODE45

MATLABodeode45

I have an ode for density that is a function of mass flow, volume, density and the time derivative of the volume. Since the volume varies as a sine function I'm expecting the densities to oscillate around the initial value, but instead it increases. Does anyone know what is causing this error and what can be done to fix it?
D0=[70 80];
tspan=[0 50];
[t T]=ode45(@DensFunc, tspan, D0);
function Dx=DensFunc(t, D) 
A=0.02;f=100;w=2*pi*f;
y=A*sin(w*t);v=w*A*cos(w*t);
a1=2e-3;
a2=1.5e-3;
VA=a1*(0.024+y);
VB=a2*(0.024-y);
dotVA=a1*v;
dotVB=-a2*v;
%setting the massflow equal to zero to find out whats wrong.  
dotmA=0;
dotmB=0;
Dx=zeros(2,1);
Dx(1)=((dotmA-D(1)*dotVA)/(VA));
Dx(2)=((dotmB-D(2)*dotVB)/(VB));
end

Best Answer

It's roundoff errors. For the large number of cycles you're integrating over, it's hard to avoid some drift; but you can lower the tolerances to get a pretty good result. Just replace
[t T]=ode45(@DensFunc, tspan, D0);
by
opts=odeset; opts=odeset(opts,'RelTol',1e-4,'AbsTol',1e-7);
[t T]=ode45(@DensFunc, tspan, D0,opts);
Related Question