Thread Subject: Solution of differential equations using ode15s

Dear Sir/Ma'am,
I am trying to solve the following multiple differential equations
using ode15s.However,I concentration of some reagents come out to be
negative.Also,the cycle stops after 5 mins,whereas it should decrease
to almost zero after about 30 mins.Please help me.

function dy=myfunc(t,y)
dy=zeros(6,1);
a=(25);
b=(10^8);
dy(1)=(13.01*y(3)*y(2)-
b*y(4)*y(1)-3*(10^7)*y(3)*y(1)-3*(10^8)*y(1)*y(2));
dy(2)=(-13.01*y(3)*y(2)-3*(10^8)*y(1)*y(2)+3.1*(10^(-3))*y(5)*y(3)+1000*y(5)*y(6));
dy(3)=(-13.01*y(3)*y(2)-
a*y(4)*y(3)-3*(10^7)*y(3)*y(1)-3.1*(10^(-3))*y(5)*y(3));
dy(4)=(-a*y(4)*y(3)-b*y(4)*y(1));
dy(5)=(13.01*y(3)*y(2)+3*(10^8)*y(1)*y(2)-3.1*(10^(-3))*y(5)*y(3)-1000*y(5)*y(6));
dy(6)=((3*(10^7)*y(3)*y(1)+3.1*(10^(-3))*y(5)*y(3)-1000*y(5)*y(6)));
end

[T,Y]=ode15s(@myfunc,[0:300:1800],[0 53900 220 3889 0 0])

Soubhik S <soubhik15@gmail.com> wrote in message <8a522366-dbcb-45e6-9f37-d5617a9769c5@a27g2000prj.googlegroups.com>...
> Dear Sir/Ma'am,
> I am trying to solve the following multiple differential equations
> using ode15s.However,I concentration of some reagents come out to be
> negative.Also,the cycle stops after 5 mins,whereas it should decrease
> to almost zero after about 30 mins.Please help me.
>
> function dy=myfunc(t,y)
> dy=zeros(6,1);
> a=(25);
> b=(10^8);
> dy(1)=(13.01*y(3)*y(2)-
> b*y(4)*y(1)-3*(10^7)*y(3)*y(1)-3*(10^8)*y(1)*y(2));
> dy(2)=(-13.01*y(3)*y(2)-3*(10^8)*y(1)*y(2)+3.1*(10^(-3))*y(5)*y(3)+1000*y(5)*y(6));
> dy(3)=(-13.01*y(3)*y(2)-
> a*y(4)*y(3)-3*(10^7)*y(3)*y(1)-3.1*(10^(-3))*y(5)*y(3));
> dy(4)=(-a*y(4)*y(3)-b*y(4)*y(1));
> dy(5)=(13.01*y(3)*y(2)+3*(10^8)*y(1)*y(2)-3.1*(10^(-3))*y(5)*y(3)-1000*y(5)*y(6));
> dy(6)=((3*(10^7)*y(3)*y(1)+3.1*(10^(-3))*y(5)*y(3)-1000*y(5)*y(6)));
> end
>
> [T,Y]=ode15s(@myfunc,[0:300:1800],[0 53900 220 3889 0 0])

You're aware that 'myfunc' doesn't use t right?

it uses t as t is the variable used to determine y...dy is actually dy/
dt

Soubhik S <soubhik15@gmail.com> wrote in message <8a522366-dbcb-45e6-9f37-d5617a9769c5@a27g2000prj.googlegroups.com>...
> Dear Sir/Ma'am,
> I am trying to solve the following multiple differential equations
> using ode15s.However,I concentration of some reagents come out to be
> negative.Also,the cycle stops after 5 mins,whereas it should decrease
> to almost zero after about 30 mins.Please help me.
>
> function dy=myfunc(t,y)
> dy=zeros(6,1);
> a=(25);
> b=(10^8);
> dy(1)=(13.01*y(3)*y(2)-
> b*y(4)*y(1)-3*(10^7)*y(3)*y(1)-3*(10^8)*y(1)*y(2));
> dy(2)=(-13.01*y(3)*y(2)-3*(10^8)*y(1)*y(2)+3.1*(10^(-3))*y(5)*y(3)+1000*y(5)*y(6));
> dy(3)=(-13.01*y(3)*y(2)-
> a*y(4)*y(3)-3*(10^7)*y(3)*y(1)-3.1*(10^(-3))*y(5)*y(3));
> dy(4)=(-a*y(4)*y(3)-b*y(4)*y(1));
> dy(5)=(13.01*y(3)*y(2)+3*(10^8)*y(1)*y(2)-3.1*(10^(-3))*y(5)*y(3)-1000*y(5)*y(6));
> dy(6)=((3*(10^7)*y(3)*y(1)+3.1*(10^(-3))*y(5)*y(3)-1000*y(5)*y(6)));
> end
>
> [T,Y]=ode15s(@myfunc,[0:300:1800],[0 53900 220 3889 0 0])

  I agree that in principle, if you start with non-negative concentrations in all six reagents, then none of them should ever become negative, as you have written the differential equations. However, with three of the reagents starting out at an actual zero concentration, which is a hair's breadth away from negative, and with some of the coefficients at very large values, it seems possible that numerical round off errors could cause some concentration, as it happens to come very close to zero, to erroneously slip into negative values. At this point it becomes possible for other concentration to also go negative.

  My suggestion would be to greatly increase the number of time points in 'tspan' from your present [0:300:1800] to say [0:10:1800] or even finer so that you can get a good plot afterwards of what happens, and then make a study of what all the concentration values are when the first one drops into negative values. That should give you an idea of what has gone wrong. You should be able to then determine whether the behavior is due to such a round off error or some other kind of malfunction.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hsmm7o$sn2$1@fred.mathworks.com>...
> I agree that in principle, if you start with non-negative concentrations in all six reagents, then none of them should ever become negative, as you have written the differential equations. However, with three of the reagents starting out at an actual zero concentration, which is a hair's breadth away from negative, and with some of the coefficients at very large values, it seems possible that numerical round off errors could cause some concentration, as it happens to come very close to zero, to erroneously slip into negative values. At this point it becomes possible for other concentration to also go negative.
>
> My suggestion would be to greatly increase the number of time points in 'tspan' from your present [0:300:1800] to say [0:10:1800] or even finer so that you can get a good plot afterwards of what happens, and then make a study of what all the concentration values are when the first one drops into negative values. That should give you an idea of what has gone wrong. You should be able to then determine whether the behavior is due to such a round off error or some other kind of malfunction.
>
> Roger Stafford

  I notice that you have not set the error tolerances in your call to 'ode15s'. That could be the source of your trouble if the concentrations are getting off to a bad start. Try setting both the relative tolerance and the six absolute tolerances using 'odeset' to appropriately tighter values than the default values call for. You have to do some experimentation with the absolute tolerances because the relative scalings among the six quantities seem greatly different. I experimented briefly with your equations and noticed that after the first few tiny steps, y(5) and y(6) both promptly became negative, using default tolerances, which theoretically they should not do.

Roger Stafford