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>
- constraint