Re: Re: [Q] Differential equation?
- To: mathgroup at smc.vnet.net
- Subject: [mg22689] Re: [mg22678] Re: [mg22620] [Q] Differential equation?
- From: Andrzej Kozlowski <andrzej at tuins.ac.jp>
- Date: Wed, 22 Mar 2000 00:27:58 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
on 3/18/00 6:27 AM, James at research at proton.csl.uiuc.edu wrote:
>>
>> James:
>>
>> These were your equations -
>>
>> y'_0(t) = -a * y_0(t) + b * y_1(t)
>> y'_1(t) = a * y_0(t) + (c*t-b) * y_1(t) --- (*)
>>
>> Interesting - what is the physical system that they describe? Can you tell
>> us?
>>
>
> Thank you for your reply.
> The system I'm working on is "2-state MMFP(Markov Modulated Fluid Process)",
> and want to compute the moment generating function of y_s(t).
> As you know, to get that function, it is necessary, as a first step,
> to solve those two equations, where I have a trouble.
> I really want to get the algebraical result, and want to know
> how I can do that.
> Thank you.
>
> James.
>
>
>
>>>
>>> Hi!
>>>
>>> I began to use Mathematica, and found out it is great.
>>> But I happen to have a question during solving differential equtations.
>>> Here's a problem.
>>>
>>> y'_0(t) = -a * y_0(t) + b * y_1(t)
>>> y'_1(t) = a * y_0(t) + (c*t-b) * y_1(t) --- (*)
>>> ^
>>> This can be solvable mathematically, even some tedious work,
>>> but when I use Mathematica, it can't solve it.
>>> After some trial and error, I found out that 't' in (*)
>>> is the problem - problem that mathematica doesn't give an answer,
>>> it just shows the above equations as an answer.
>>> So I wonder if this is the limit of Mathematica,
>>> or is there any way to solve it?
>>> I sincerely hope there's some way - because my work involves
>>> a lot of Diffrential Equations.
>>> Any reply would be appreciated.
>>>
>>>
>>> James.
>>>
>>>
>>
>>
>>
>
>
>
A complete algebraic solution has already been sent by Jens-Peer Kuska.
However, it is is my impression that it was not written in a way that would
be easy to follow for someone relatively new to Mathematica and so I thought
it might be worthwhile to try to give a more detailed presentation of it.
To start with, it is not really surprising that Mathematica does not solve
this system of equations automatically because it is a system of ordinary
differential equations with non-constant coefficients. In the case of
constant coefficients there is a perfectly general algorithm but when the
coefficients are not constant only certain special cases can be solved in
terms of explicit formulas. In most of such cases human intervention at some
point seems to be necessary.
Let us start by defining two expressions eq1 and eq2 such as our equations
become eq1==0 and eq2==0;
In[1]:=
eq1 = y[1]'[t] - a*y[0][t] - (c*t - b)*y[1][t];
In[2]:=
eq2 = y[0]'[t] + a * y[0][t] - b * y[1][t];
We shall add to these a third one obtained by differentiatig one of them,
lets say eq1:
In[3]:=
eq3 = D[eq1, t]
Out[3]=
-(c y[1][t]) - a (y[0])'[t] - (-b + c t) (y[1])'[t] + (y[1])''[t]
We can now elliminate y[0],and y'[0] from the equations to obtain a single
equation of second order with non-constant coefficients. One can do this
elimination by hand but it is more convenient to use GroebnerBasis:
In[4]:=
gr1 = GroebnerBasis[{eq1, eq2, eq3}, {y[1]''[t], y[1]'[t],
y[1][t]}, {y[0]'[t], y[0][t]}, MonomialOrder -> EliminationOrder]
Out[4]=
{-(c y[1][t]) - a c t y[1][t] + a (y[1])'[t] + b (y[1])'[t] - c t (y[1])'[t]
+ (y[1])''[t]}
We can now get our second order differential equation as follows:
In[5]:=
diff1 = Collect[gr1[[1]], {y[1]''[t], y[1]'[t], y[1][t]}]
Out[5]=
(-c - a c t) y[1][t] + (a + b - c t) (y[1])'[t] + (y[1])''[t]
Mathematica is unable to solve this equation with
DSolve[diff1==0,y[1][t],t]. However, there is a well known method of
converting any linear differential equation of the form
y''[t]+P[t]*y'[t]+Q[t]*y==0 to the "standard form" z''[t]+q[t]*z[t]==0.
Equations in this standard form have a better chance of being solved (though
even this is not guaranteed). Here is one way to obtain an equation in
standard form:
In[6]:=
diff2 = diff1 /. y[1] -> (E^(Log[z[#]] + (1/2)Integrate[-(a + b - c #), #])
&) // Simplify
Out[6]= -((E^((t*(-2*a - 2*b + c*t))/4)*((a^2 + b^2 - 2*b*c*t + 2*a*(b +
c*t) + c*(2 + c*t^2))*z[t] - 4*Derivative[2][z][t]))/ 4)
Observe that there is no z'[t] term. Mathematica really now ought to be able
to manage:
DSolve[diff2==0,z[t],t]
but it still fails to produce an answer. This is strange because it turns
out that this is due to the E^something factor in the expression diff2 which
of course is quite irrelevant to the equation diff2==0. Once we get rid of
it Mathematica will at last produce a solution:
In[7]:=
diff3 = diff2 /. E^_ -> -4
Out[7]=
(a^2 + b^2 - 2*b*c*t + 2*a*(b + c*t) + c*(2 + c*t^2))*z[t] -
4*Derivative[2][z][t]
This is the same equation as in Jens Kuska's message.
In[8]:=
sol1 = DSolve[diff3 == 0, z[t], t] // FullSimplify
Out[8]=
{{z[t] -> (((a - b + c*t)^2)^(3/4)*(C[2]*Hypergeometric1F1[(3 +
Sqrt[c^(-2)]*(2*a*b + c))/4, 3/2,
(Sqrt[c^(-2)]*(a - b + c*t)^2)/2] + C[1]*HypergeometricU[(3 +
Sqrt[c^(-2)]*(2*a*b + c))/4, 3/2,
(Sqrt[c^(-2)]*(a - b + c*t)^2)/2]))/(2*E^((Sqrt[c^(-2)]*(a - b +
c*t)^2)/4)*(-(c^2*(a - b + c*t)^2))^(1/4))}}
The last step is to express to substitute back to get a solution of the
original equation:
In[9]:=
(E^(Log[z[t]] + (1/2)Integrate[-(a + b - c t), t])) /. sol1[[1, 1]]
Out[9]=
E^(((-a - b)*t + (c*t^2)/2)/2 +
Log[(((a - b + c*t)^2)^(3/4)*(C[2]*Hypergeometric1F1[(3 +
Sqrt[c^(-2)]*(2*a*b + c))/4, 3/2,
(Sqrt[c^(-2)]*(a - b + c*t)^2)/2] + C[1]*HypergeometricU[(3 +
Sqrt[c^(-2)]*(2*a*b + c))/4, 3/2,
(Sqrt[c^(-2)]*(a - b + c*t)^2)/2]))/(2*E^((Sqrt[c^(-2)]*(a - b +
c*t)^2)/4)*(-(c^2*(a - b + c*t)^2))^(1/4))])
This "ought to" be the general solution. I say "ought to", because I would
be very careful when using it for any practical purpose. The function
HypergeometricU has a branch cut along the negative real axis. While I am
not an expert in this area to me this suggests that the solution may not be
correct for all values of the parameters. One should investigate this by
substituting numerical values anc checking if the functions thus obtained
satisy the original equation.
--
Andrzej Kozlowski
Toyama International University
JAPAN
http://sigma.tuins.ac.jp