Preface

This is a tutorial made solely for the purpose of education and it was designed for students taking Applied Math 0330. It is primarily for students who have very little experience or have never used Mathematica before and would like to learn more of the basics for this computer algebra system. As a friendly reminder, don't forget to clear variables in use and/or the kernel.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in normal font. This means that you can copy and paste all commands into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

Numerical Solutions

Sometimes Mathematica is able to solve a differential equation. However, it is a rear case and in practice, we rely on numberical methods. Most differential equations that occur in actual practice are much too complicated to admit a solution expressed through elementary or special functions. Usually, we are happy to obtain a solution in implicit form, but unfortunately it happen ocationally. The fact that a differential equation does not possess an explicit or implicit solution does not mean that the solution does not exist.

Since the kernel of a differential operator usually contains a space of functions, we do not expect to obtain a general solution that includes arbitrary constants. Instead, we knock out one solution by considering the initial value problem (IVP for short)

Example: Consider the initial value problem \( y' = \sin \left( x^2 \right) , \quad y(0) =0 . \) We try to find its solution using standard Mathemaica command DSolve:

exactsol = DSolve[{y'[x] == Sin[x^2], y[0] == 0}, y[x], x]

Out[1]= {{y[x] -> Sqrt[\[Pi]/2] FresnelS[Sqrt[2/\[Pi]] x]}}

Mathematica represents the solution through the special function, called Fresnel function: \( S(x) = \int_0^x \sin \left( t^2 \right) {\text d}t . \)

?FresnelS

FresnelS[z] gives the Fresnel integral S(z).  >>

Now we solve the same problem numercally using NDSolve.

numsol = NDSolve[{y'[x] == Sin[x^2], y[0] == 0}, y[x], {x, 0, 10}]

Out[3]= {{y[x] -> InterpolatingFunction[{{0., 10.}}],<>],[x]}}

Plot[y[x] /. numsol, {x, 0, 10}, PlotRange -> All]

This is a plot of the numerical solution. Now we calculate its approximate value at x=6.5 as follows:

numsol /. x -> 6.5

Out[5]= {{y[6.5] -> 0.639918}}

This indicates that y(6.5) approximately equal to 0.639918. We can improve accuracy by choosing the option PrecisionGoal

There are two main approaches to find a numerical value for the solution to the initial value problem \( y' = f(x,y), \quad y( x_0 )= y_0 \) when its solution can be obtained using either DSolve (for solution represented via known functions, if it is possible) or NDSolve (for numerical solution) command. The first one is based on definition of the slope function f(x,y) as an expression. The second one (which is recommended) is to define the slope as a function. To see how it works, let us consider an example.

 As we see application of NDSolve is very straightforward --- use the solution "interpolating function" as a regular function. When plotting solution, Mathematica recommends to use Evaluation command. We give some examples.

Example: : First, we solve a first-order ordinary differential equation, which gives you sort of decaying oscillations:

sol1 = NDSolve[{y'[x] == y[x] Sin[2 x + y[x]], y[0] == 1}, y, {x, 0, 30}];

Now we use it inside another equation. The important thing to remember is that the domain of the previous equations should be equal or contain inside the domain of next equations.

sol2 = NDSolve[{z'[x]== -z[x]^2+First[Evaluate[y[x] /. sol1]], z[0] == 0}, z, {x, 0, 30}]

Then we plot both solutions. Indeed, first solution here plays the role of sort of a driving force, so on the graphs below you see synchronization between oscillations.

Plot[{Evaluate[y[x] /. sol1], Evaluate[z[x] /. sol2]}, {x, 0, 30}, PlotRange -> All]

Now you can do more advanced things, like using also derivative of the "Interpolating Function" inside of new equations:

sol3 = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}];
sol4 = NDSolve[{z'[x]== -z[x]+First[y[x]/5 + y'[x] /. sol3], z[0] == 0}, z, {x, 0, 30}];
Plot[{Evaluate[y[x] /. sol3], Evaluate[z[x] /. sol4]}, {x, 0, 30}, PlotRange -> All, Filling -> 0]