Re: constraint

*To*: mathgroup at smc.vnet.net*Subject*: [mg32758] Re: [mg32741] constraint*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Fri, 8 Feb 2002 03:49:33 -0500 (EST)*References*: <200202071012.FAA07513@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

sahni mayank wrote: > > i'm having problems with inserting an algebraic > constraint into a series of differential equations and > solving with NDSolve. > for example > sol1 = NDSolve[{ > y1'[t] == 0.5*k2 + k4*y2[t]*y3[t] - > k10*y1[t]*y4[t], > y2'[t] == > k1 + k3 + k5*y1[t]*y1[t] - k4*y2[t]*y3[t], > y3'[t] == > k9*y3[t]*y1[t] - k4*y2[t]*y3[t] - > k8*y1[t]*y3[t], > > and instead of using a differential eq for y4[t] i > want to use an algebraic equation for example > y4[t]=k12-y1[t]-y2[t], > > is it possible to solve this using NDSolve? > > __________________________________________________ > Do You Yahoo!? > Send FREE Valentine eCards with Yahoo! Greetings! > http://greetings.yahoo.com This topic is discussed at the link below. http://support.wolfram.com/mathematica/mathematics/equations/dae.html One method is to eliminate the algebraic equations(s). We do this for your example, with numeric coefficients set at random. rul = Thread[{k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12}-> Table[Random[],{12}]]; sol1 = NDSolve[{ y1'[t] == 0.5*k2 + k4*y2[t]*y3[t] - k10*y1[t]*y4[t], y2'[t] == k1 + k3 + k5*y1[t]*y1[t] - k4*y2[t]*y3[t], y3'[t] == k9*y3[t]*y1[t] - k4*y2[t]*y3[t] - k8*y1[t]*y3[t], y1[0]==1, y2[0]==3, y3[0]==-2} /. y4[t]->(k12-y1[t]-y2[t]) /. rul, {y1[t],y2[t],y3[t]}, {t,0,1}] When such elimination is difficult, one might instead convert to a pure ODE system. This is easy to do just by differentiating the algebraic constraint(s). Below I show this for your example. Note that initial conditions must satisfy relevant algebraic constraints. sol2 = NDSolve[{ y1'[t] == 0.5*k2 + k4*y2[t]*y3[t] - k10*y1[t]*y4[t], y2'[t] == k1 + k3 + k5*y1[t]*y1[t] - k4*y2[t]*y3[t], y3'[t] == k9*y3[t]*y1[t] - k4*y2[t]*y3[t] - k8*y1[t]*y3[t], y4'[t] == -y1'[t]-y2'[t], y1[0]==1, y2[0]==3, y3[0]==-2, y4[0]==k12-y1[0]-y2[0]} /. rul, {y1[t],y2[t],y3[t],y4[t]}, {t,0,1}] You can assess quality of the result as below. Plot[(y4[t] - (k12-y1[t]-y2[t]) /. rul) /. sol2[[1]], {t,0,1}] You can also do something like Plot[(y1[t]/.sol1) - (y1[t]/.sol2), {t,0,1}] to see that the solutions are in agreement. One can often improve on this second method in cases where the solution wanders off the algebraic constraint variety. The basic idea is to chop up the time interval and use a root-finder such as FindRoot to move back to that variety as a way to correct the missteps. Daniel Lichtblau Wolfram Research

**References**:**constraint***From:*sahni mayank <sahnimayank@yahoo.com>