How can I solve a system of differential equations with a specific form?

[cesca]

Hi,

I am a chemist with unfortunately too little experience with differential equations, and now I have a problem I need to solve:
I have a system of DE of the form:

y' = ay-b*y^2-c*y*x
x'=d*y^2

a, b, c, and d are constants (related to each other)

I have been trying to read about this, but I got stuck after having read about first order systems. I would very much appreciate any help, and I hope that I am posting in the right forum.

 

[Odysseus]

Hi cesca,

We'll need a little more information: initial conditions, the relationships between a, b, c and d, the interval where you hope to find solutions. Out of curiosity, I'd be interested in knowing the system you are attempting to model. I'm assuming all derivatives are with respect to t?

 

[JJacquelin]

y' = ay-b*y^2-c*y*x
x'=d*y^2
supposing that y = y(t)

You can find a particular solution on the form :
y = k/t
x = h -d*k²/t
k and h constants to be computed (related to a, b, c, d)

 

[jackmell]

Hi,

I am a chemist with unfortunately too little experience with differential equations, and now I have a problem I need to solve:
I have a system of DE of the form:

y' = ay-b*y^2-c*y*x
x'=d*y^2

a, b, c, and d are constants (related to each other)

I have been trying to read about this, but I got stuck after having read about first order systems. I would very much appreciate any help, and I hope that I am posting in the right forum.

Hi. Nice if we explicitly state "analytic" solution or just anyway of solving it. You didn't so you can do a very nice job of solving it numerically if that is acceptable. Just design an initial value problem with all the constants assigned numeric values and y(0) and x(0) and then just run it through the DE numeric solver in Mathematica: NDSolve. In fact, you could do that, then fit x(t) and y(t) to some polynomials and then those polynomial approximations would be a nice "fit" to the analytic solution again, if that is acceptable to you:

a = 1; 
b = 1; 
c = 1; 
d = 1; 
x0 = 1; 
y0 = 2; 
mysol = NDSolve[{Derivative[1][y][t] == 
     a*y[t] - b*y[t]^2 - c*y[t]*x[t], 
    Derivative[1][x][t] == d*y[t]^2, 
    y[0] == y0, x[0] == x0}, {x, y}, 
   {t, 0, 1}]
myplot = Plot[{x[t], y[t]} /. mysol, 
   {t, 0, 1}]
myxdata = Table[
   {t, First[x[t] /. mysol]}, 
   {t, 0, 1, 0.1}]
myxfit = Fit[myxdata, Table[t^n, 
    {n, 0, 10}], t]
myydata = Table[
   {t, First[y[t] /. mysol]}, 
   {t, 0, 1, 0.1}]
myyfit = Fit[myydata, Table[t^n, 
    {n, 0, 10}], t]
myfitplot = Plot[{myxfit, myyfit}, 
   {t, 0, 1}]
Show[{myplot, myfitplot}]

And if the system is well-behaved, the 10-degree polynomial functions myxfit[t] and myyfit[t] are probably now very good "analytic" approximations to the solution for this particular IVP in the range [0,1].

 

[JJacquelin]

Sure, Jackmell is right : numerical calculus is tne best method on a practical viewpoint.
Nevertheless, analytical solution is possible (attachment). The result is obtained on a parametric form but the formulas are rather complicated.

 

[jackmell]

Sure, Jackmell is right : numerical calculus is tne best method on a practical viewpoint.
Nevertheless, analytical solution is possible (attachment). The result is obtained on a parametric form but the formulas are rather complicated.

Hi Jacquelin. I'm going through your solution so I can compare it to the numerical solution.

You have:

Let:

$$X=Cx-A;\quad dX=Cdx\rightarrow DC\frac{Dy}{dX}=\frac{X}{y}-B$$

should that not be:

$$X=Cx-A;\quad dX=Cdx\rightarrow DC\frac{Dy}{dX}=-\frac{X}{y}-B$$

 

[JJacquelin]

You are right, that's a careless mistake. May be not the only one !
The point is the method of resolution.
Well, to be corrected and more important, to be verified.

 

[cesca]

Wow - a lot of answers. I will try to go trhough them and ask again if I get stuck (which I probably will). A numerical solution is perfectly acceptable, expecially if the analytical one is very complicated (or not possible to get).

Odysseus is right, it is a time-derivative, and I am looking at how the concentration of some chemical species change over time. I have a suggestino for a reaction mechanism, and sould like to know if that suggestion fits with my experimental data.
Relationships between a, b, c, and d are:

a = y*k*OHs
b = y*(k+l)
c = y*k
d= y*l

Initial conditions: I am not entirely sure, [hc], should be small but not zero, since this reaction mechanism would not be the predominant one before some [hc] has been buildt up.

Thank a lot!

 

[jackmell]

Hi. I've compared Jacquelin's method against the numerical results although I avoided a value that would give X=0 since we're taking logs below. So I let:
$$a=1$$
$$b=1$$
$$c=2$$
$$d=1$$
$$x_0=1$$
$$y_0=2$$

I believe there is a mistake in the last line of Jac's analysis. Starting with the expression:

$$\frac{dX}{X}=-\frac{DCYdY}{1+BY+CDY^2}$$

then:

$$\ln(x)\biggr|_{X_0}^{X}=J(Y)-J(Y_0)$$

where J(Y) is the antiderivative on the right side and [itex]X_0=Cx_0-a[/tex] and [itex]Y_0=\frac{y_0}{X_0}[/tex]

X(Y) and Y are as before but I think t is calculated this way unless I'm missing something:

$$t=\frac{1}{D}\int \frac{dx}{y^2}=\frac{1}{CD}\int \frac{d(X(Y))}{X(Y)^2 Y^2}$$

This gives a positive value of t for Y in the range of (0,2].

Here's the code I used to compare the numerical results to the analytic results. The final plot superimposes the two methods and they agree very well in the range I computed them although there are some further considerations for Y very close to zero.

a = 1; 
b = 1; 
c = 2; 
d = 1; 
x0 = 1; 
y0 = 2; 
mysol = NDSolve[{Derivative[1][y][t] == a*y[t] - b*y[t]^2 - c*y[t]*x[t], 
    Derivative[1][x][t] == d*y[t]^2, y[0] == y0, x[0] == x0}, {x, y}, 
   {t, 0, 5}]
myplot = Plot[{x[t], y[t]} /. mysol, {t, 0, 5}]
myantid[Y_] := (1/2)*((2*b*ArcTan[(b + 2*c*d*Y)/Sqrt[-b^2 + 4*c*d]])/
     Sqrt[-b^2 + 4*c*d] - Log[1 + b*Y + c*d*Y^2])
Y0 = 2; 
theX[Y_] := Exp[myantid[Y] - myantid[Y0]]
mydx = FullSimplify[D[theX[u], u]]
myx[Y_] := (theX[Y] + a)/c; 
myy[Y_] := theX[Y]*Y; 
thet[Y_] := (1/(c*d))*NIntegrate[mydx/(theX[u]^2*u^2), {u, Y0, Y}]
pp1 = ParametricPlot[{{thet[Y], myx[Y]}, {thet[Y], myy[Y]}}, 
   {Y, 0.001, 2}, PlotPoints -> 100]
Show[{myplot, pp1}]

This was a fun problem and I think a very good example of manipulating differentials. Thanks Jacquelin. :)