Today, I tried some more real world examples. May plan is to solve the fall of a body with air resistance according to

\(y“ = -g+k \, y’^2\)

using g=9.81, k=0.01 (with IS units), and a drop height of 100. The question to answer is: When does the body hit the ground?

To learn about the Mathematica command DSolve, I tried

DSolve[y'[x]=x+y[x],y[x],x]

which does not work. The reason is that an equation writes with == not =. I noticed that, and tried

DSolve[y'[x]==x+y[x],y[x],x]

which also does not work! The reason is now that y'[x] has been assigned a value by my error above. This is really tough, I must say!

To recover, I tried „Clear[y]“ and „Clear[y‘], but that does not work. However, „y'[x]:=y'[x]“ seemed to work. I then learned, how to extract the function from the solution, and how to replace the constant with some value.

solution  = DSolve[y'[x] == x + y[x], y[x], x]
f = y[x] /. solution
Plot[f /. C[1] -> 1, {x, -1, 1}]

This produces the correct plot.

Back to our problem. It worked out nicely in Mathematica. The initial problem gets an immediate solution. In fact, three solutions, two of which involve the complex unit i.

We need the first solution. So we extract it and plot it with time x from 0 to 8.

f = y[x] /. solution[[1, 1]]
Plot[f /. {g -> 9.81, k -> 0.01, h -> 100}, {x, 0, 8}]

It is now easy to solve this for  the value 0 (unless you forget to use == and spoil everything with =).

Solve[(f /. {g -> 9.81, k -> 0.01, h -> 100}) == 0, x]
  {{x -> -5.29184}, {x -> 5.29184}}

Not all differential equations can be solved analytically. E.g., if you enter

DSolve[y''[x] == -g + k*Abs[y'[x]]^1.67, y[x], x]

you get no answer at all. In fact, Mathematica will run and run, and run. To stop, press Alt+Dot. What we need is a numerical solution. The usual method, which I found in a book about Mathematica, returns an interpolating function.

Unfortunately, I could not find a way to solve this for 0. The answer is in terms of InverseFunction, which is no help at all.

Any numerical solver of a program like Euler or Matlab should not have troubles with this equation. Let us try. We have to rewrite this into a system of equations.

>function f(x,y) := [y[2];-g+k*abs(y[2])^a]
>g=9.81; k=0.01; a=1.67;
>x=0:0.01:8; y=ode("f",x,[100;0]);
>plot2d(x,y[1]):

To solve for 0, we can use the numerical output and splines, or write a function, which solve the equation, and use the secant method on this function.

>s=spline(x,y[1]);
>function fsol(t) := evalspline(t,x,y[1],s)
>solve("fsol",4.5)
  4.78765346167
>function g(t) := ode("f",[0,t],[100;0])[1,2]
>solve("g",4.5)
  4.78765346143

How about the exact solution of the problem with exponent 2 in Maxima? It turns out, that Maxima can find the general solution, but it cannot solve for the constants. I did not find an obvious trick to overcome this.

To find a solution of the initial value problem, you have to solve for y‘ first. The solution for y‘ is implicit, and needs to be solved first, which is no problem. You can solve for the constant so that the initial condition y'(0)=0 is met. Then you integrate y‘.

Here are the command with no printed output.

>&assume(g>0,k>0);
>sol &= ode2('diff(yd,x)=-g+k*yd^2,yd,x);
>ydsol &= rhs(solve(sol,yd)[1]);
>&solve((ydsol with x=0)=0,%c); ydsol &= ratsimp(ydsol with %[1]);
>&assume(t>0); function f(t) &= integrate(ydsol,x,0,t)+h;
>g=9.81; k=0.01; h=100;
>solve("f",5)
  5.29184480915

This, of course, is not a nice solution.

What about Maple? My first attempt returned a solution involving RootOf, even though I told Maple to assume g>0 and k>0. This is correct, and you can get rid of it with allvalues. Here is a screen dump.

Then

solve(g=0,x)

returns the solution. The numerical solution can be obtained with the type=numeric flag. It returns a procedure, which can be plotted with the function odeplot. It can also be evaluated at special points. But I found no way to solve it for 0.