Using dsolve instead of solve

JacquesC's picture

solve used to be one of Maple's strongest commands -- it even subsumed simplify in power.  But, over the years, dsolve slowly took over as the most powerful comand.  At the same time, people started realizing that within the framework of differential equations, the toolbox was actually larger than the one for algebraic equations (and most algebraic tools are still available).  So many tasks that one thinks of doing purely algebraically can also be done using differential equations, with perhaps the most surprising one is to factor multivariate polynomials via partial differential equations.

Here we show how to use dsolve to solve parametric equations, ie equations that depend on a parameter -- the availability of a parameter is crucial for this method.  As a first example, let us look at solving a simple version of a cubic equation: x^3+b*x+a=0.  This actually has 2 parameters, so let us choose a as the one to concentrate on.  We want to express this equation as a differential equation -- and luckily there is a command for just that, namely gfun[algeqtodiffeq].  We can pass the result directly to dsolve.

eq1 := x^3 + a*x + b = 0:
de := gfun[algeqtodiffeq](eq1,x(b)):
dsolve(de);

and we get the result as

                       2    3
  x(b) = -1/3 RootOf(_Z  + a , index = 1)

                        2      2         1/2
        (3 RootOf(_Z (_Z  + a))  + 2 a) 3

                               1/2
                            3 3    b          /  2
        sin(1/3 arctan(------------------))  /  a
                            3       2 1/2   /
                       (-4 a  - 27 b )

                                                       1/2
                        2                           3 3    b
         + RootOf(_Z (_Z  + a)) cos(1/3 arctan(------------------))
                                                    3       2 1/2
                                               (-4 a  - 27 b )

The first thing to notice is that this is not Cardano's formula, but the much-better trigonometric version of the cubic formula! The second thing we notice is that we can ``see'' the 3 different branches , and the choice that depends on the initial value of a.

Why show this technique now, other than it is fun?  Because of an older post on solve, asking to solve a exp+trig problem.  We can do the same thing here, but we have to use PDEtools[dpolyform] instead.  The command sequence is then

eq3 := exp(-1/k*t)*k/(1+k^2)+(-cos(t)*k+sin(t))/(1+k^2):
ee := eval(eq3,t=t(k)):
de := PDEtools[dpolyform](ee=0,t(k),no_Fn):
dsolve(op([1,1],de));

which does give a closed-form, with 2 arbitrary constants (as the DE was of order 2). We need initial conditions to get rid of those constants. We can work out easily enough that t(0)=0, and with a little more work that D(t)(0) = RootOf(1-exp(_Z)+_Z*exp(_Z)), which is in fact LambertW(_NN1, -exp(-1))+1. Passing that, as is, as initial conditions to dsolve, I gave up after a few minutes. But with some perseverance, it is possible that one might be able to find a closed-form for this equation!

 

Comments

jpmay's picture

Factoring via PDEs

Shuhong's paper on multivariate factoring via PDEs is great, but a slighly funny example here.  The PDE in that paper is apparently pretty strange (notwithstanding the fact that it is considered as being defined over finite fields!), and for the purposes of the algorithm it is instantly turned into a matrix - so it is not actually using any of the theory of PDEs.  I suspect most of the reason for putting 'PDE' in the title is was to make it eye catching.

John

JacquesC's picture

You're right

I should have found a better example (as there are so many).  The title of the paper was simply too good for me to pass up...

Robert Israel's picture

Not very closed form

What dsolve gives me (in Maple 11) is in no way "closed form": it involves a big RootOf an expression involving logarithms and trig functions, and simplifies (with combine) to

t(k) = ln(k)*k-1/2*ln(1+k^2)*k+1/2*k*ln(1/RootOf(k*sin(-1/2*ln(1/_Z^2)*k+_C2*k
-_C1-ln(k)*k+1/2*ln(1+k^2)*k)-cos(-1/2*ln(1/_Z^2)*k+_C2*k-_C1-ln(k)*k+1/2*ln(1
+k^2)*k)+_Z*(1+k^2)^(1/2))^2)-_C2*k

 I don't see this as an improvement over the result of solve, even if the values of _C1 and _C2 could be found.  By the way, k=0 might not be a good place to take the initial conditions, because of the singularity there in the differential equations (and presumably in the solution).

jpmay's picture

Solve doing all the work?

This is an interesting idea.

However, in this case, it is entirely possible that solve is even doing "all" the work.  Putting a trace on solve before executing these commands shows, for example, solve is being called:

solve(ln(x)-1/2*ln(1+x^2)-1/x*y(x)+
   1/2*ln(1/cos(y(x)-arctan(x)+_C1)^2)-_C2, {y(x)});
solve( x/(1+x^2)^(1/2)/exp(1/x*y(x))*(1/(cos(y(x))/(1+x^2)^(1/2)
   *cos(_C1)+cos(y(x))*x/(1+x^2)^(1/2)*sin(_C1)+sin(y(x))
   *x/(1+x^2)^(1/2)*cos(_C1)-sin(y(x))/(1+x^2)^(1/2)
   *sin(_C1))^2)^(1/2)/exp(_C2)-1, y(x));

And those seem almost as complicated as the original problem.

John

JacquesC's picture

In this case

In the problem that reminded me of this method, you are probably right that solve is in fact doing all the hard work.  Which is good, in a way - it means that dsolve is properly calling solve!

The fundamental technique works rather well though, in many cases yielding solutions that are interestingly different than solve's.  To me, this says that something 'deeper' is going on -- that on some problem, the 'smooth' structure contains information which can somehow be leveraged in different ways than the 'algebraic' structure. 

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
}