Obtaining the Roots of Non-Linear Equations and Integrating Some Functions

 

This article is designed for those interested in the advanced topics in Mathematica, dealing with solving complex equations and graphing those equations. The Built-in- Function “FindRoot” and “Plot” are used to solve complex equations associated with Richard’s equations and graphing them. These equations appear in the field of Soil Physics and in other math and engineering applications. In the following, some equations involving Error Function and the Complementary Error Function are graphed using Built- in-Function “ FindRoot”. Also, the integral of a few famous functions are obtained and the results are compared with those in the integral tables.

 

 

Use the routine “FindRoot” to find the Roots of the following equation:

 

FindRoot[2.18*((t/p)^.5)*Exp[-t/4.] -.216==0,{t,.005}]

 

And the answer is:

 

{0®0.0313291}

 

However, if we change the search for the roots from .005 to 20, we arrive at a new root, here is how:

 

 

FindRoot[2.18*((t/p)^.5)*Exp[-t/4.] -.216==0,{t,20}]

 

 

And the answer is:

 

{0®11.913}

 

Built-in-Function “FindRoot” is an excellent means for finding the roots of non-linear equations.

 

Users may notice that in order to solve a non-linear equation such as this we should have some ideas about the location of the desired root – in this case “11.913”. This is equation (13) from the paper by Chen et al. published in Water Resources Research vol. 37, No. 4. pages 1091-1093. Equation (13) represents surface flux after the soil parameters are substituted into equation (13). Now we are going to plot the above equation to get some ideas about its behavior, like this:

 

Plot[2.18*((t/p)^.5)*Exp[-t/4.] -.216,{t,0,15}]

 

and the graph produced using above plot procedure is:

 

 

 

 

 

To better observe the behavior of the root near zero, t=0.03132

 

 

 

We change the limits to:

 

Plot[2.18*((t/p)^.5)*Exp[-t/4.] -.216,{t,0,1}]

 

 

And the root near zero now appears in the above graph:

 

In the above graph, it appears that the curve crosses the x axis somewhere below (.05), we observed earlier that the exact solution is: t =.031329.

 

In the following equation t, time is set to t=100 second and the equation is graphed for z.

 

Plot[ .45 -.00667*Exp[(z/2)-.00000055]*(((Exp[-z^2/.0000088])/1.77)-(z/.00297)*Erfc[z/.000297]),{z,0.0,1.6}]

 

 

 

 

 

 

Below, the numerical value for x-axis is changed from 1.6 to .2, so the plot may cover a smaller but more accurate domain and the above graph looks clearer.

 

 

In the following equation time, t is equal to two million seconds and is another plot example of the complementary Error Function.

 

Plot[  .45   -.94 * Exp [ ( t /2 ) -.11] * ( ( ( Exp [ - t ^2 /.18] ) /1.77 ) - ( t /.42 ) * Erfc [ t /.42] ) , {t,0.0,1.6} ]

 

 

 

 

 

 

 

 

Below time, t assumes a small number and plot routine is designed to yield two curves on the same axis.

 

tt={.45   -.00667 * Exp [ ( t /2 ) -.00000055] * ( ( ( Exp [ - t ^2 /.0000088] ) /1.77 ) - ( t /.00297 ) * Erfc [ t /.000297] ) ,.45   -.94 * Exp [ ( t /2 ) -.11] * ( ( ( Exp [ - t ^2 /.18] ) /1.77 ) - ( t /.42 ) * Erfc [ t /.42] )   } ;

 

pp=Plot[Evaluate[tt],{t,.0,1.6}];

 

resulting in:

 

 

 

 

 

 

 

 

In the following, the Built-in-Function “Plot” includes five similar but slightly different equations, which results in five graphs on the same x-y axis.

 

tt={.45       -.00667     *     Exp     [     (     t     /2     )     -.00000055]     *     (     (     (     Exp     [     -     t     ^2     /.0000088]     )     /1.77     )     -     (     t     /.00297     )     *     Erfc     [     t     /.000297]     )     ,.45       -.211     *     Exp     [     (     t     /2     )     -.00055]     *      (   (     (     Exp     [     -     t     ^2     /.0088]     )     /1.77     )     -      (    t     /.09     )     *   Erfc     [     t     /.09]     )   ,

      .45       -.472     *     Exp     [     (     t     /2     )     -.0027]     *     (     (     (     Exp     [     -     t     ^2     /.044]     )     /1.77     )     -     (     t     /.209     )     *     Erfc     [     t     /.209]     )   ,.45       -.670     *     Exp     [     (     t     /2     )     -.0055]     *     (     (     (     Exp     [     -     t     ^2     /.088]     )     /1.77     )     -     (     t     /.30     )     *     Erfc     [     t     /.30]     )   ,

      .45       -.94     *     Exp     [     (     t     /2     )     -.011]     *     (     (     (     Exp     [     -     t     ^2     /.18]     )     /1.77     )     -     (     t     /.42     )     *     Erfc     [     t     /.42]     )       }     ;

 

pp=Plot[Evaluate[tt],{t,0.0,1.6}];

 

 

and the graphs is :

 

 

 

 

Please note that in the above five graphs are created by Plot, the one with a smaller (close to zero) x-value almost rests on y-axis. Also note that the fourth graph from the left goes over the line x=.5 this is not correct and is due to a small error made. In the above equation if .3 is changed to .2966 (which is the correct value) then the result will not cross the x=.5 line. This shows the sensitivity of the Error Function:

 

In the following Plot the correction is made to the above equation:

 

tt={.45         -.00667       *       Exp       [       (       t       /2       )       -.00000055]       *       (       (       (       Exp       [       -       t       ^2       /.0000088]       )       /1.77       )       -       (       t       /.00297       )       *       Erfc       [       t       /.000297]       )       ,.45         -.211       *       Exp       [       (       t       /2       )       -.00055]       *        (     (       (       Exp       [       -       t       ^2       /.0088]       )       /1.77       )       -        (      t       /.09       )       *     Erfc       [       t       /.09]       )     ,  

       .45         -.472       *       Exp       [       (       t       /2       )       -.0027]       *       (       (       (       Exp       [       -       t       ^2       /.044]       )       /1.77       )       -       (       t       /.209       )       *       Erfc       [       t       /.209]       )     ,.45         -.670       *       Exp       [       (       t       /2       )       -.0055]       *       (       (       (       Exp       [       -       t       ^2       /.088]       )       /1.77       )       -       (       t       /.2966       )       *       Erfc       [       t       /.2966]       )     ,  

       .45         -.94       *       Exp       [       (       t       /2       )       -.011]       *       (       (       (       Exp       [       -       t       ^2       /.18]       )       /1.77       )       -       (       t       /.42       )       *       Erfc       [       t       /.42]       )         }       ; 

 

pp=Plot[Evaluate[tt],{t,0.0,1.6}];

 

 

 

 

 

 

To graph the integral of the Complementary Error Function below, from 0.0 to 1.6:

 

 

The left hand side of the above equation is the integral of the Error Function.

First, the dimensionless Time, T, is set to T =1/4, and after putting the above equation into the Mathematica’s format we arrive at:

 

 

tt=.564Exp[-z^2]-z*Erfc[z];

 

pp=Plot[Evaluate[tt],{z,0.0, 1.6}]

 

And the graph of the integral of the Complementary Error Function is:

 

 

 

 

 

Now graph the integral of the above integral or graph the function below:

 

 

Again after letting T=1/4, simplifying the above equation and putting it in the Mathematica’s format we arrive at:

 

tt=-.282*z*Exp[-z^2]+((1.+2z^2)/4)*Erfc[z];

 

pp=Plot[Evaluate[tt],{z,0.0, 1.6}]

 

And the graph of double integral of the Complementary Error Function is: 

 

 

The above two equations for the two integral of the Complementary Error Function are placed under one equation and one graph. Below is the equation:

 

 

tt={.564Exp[-z^2]-z*Erfc[z],-.282*z*Exp[-z^2]+((1.+2z^2)/4)*Erfc[z]};

 

pp=Plot[Evaluate[tt],{z,0.0, 1.6}]

 

and the answer is:

 

 

 

 

It is left to the readers to name the axis and identify the graphs, which one belongs to iErfc(z), and which one to i^2* of Erfc(z).

 

Now lets integrate the following:

 

 

 

and the result is:

 

 

 

To obtain the integral of:

 

 

 

Please note that the above integral is very similar to the integral number 3.6214 on page 369 of the integral book by Gradshteyn and Ryzhik. 

 

 

 

Now lets look at some fancy integration using Mathematica:

 

Integrate the following integral:

 

Integrate[(Cos[x])^(2m + 1), {x, 0, p/2}]

 

And the result is:

 

In the above integral if “p/2” is changed to “π/2”, then the result is:

 

 

Please note that the above is the same integral given in “Table of Integrals ….” By Gradshteyn & Ryzhik, known as the Russian Integral Book.

 

 

We conclude our examination of Mathematica by solving and graphing equation [114] on pages 835 of the paper by Ron Wiltshire et al. published in Journal of Physics. This paper is one of the best papers published in the filed of Soil Physics in the J. Phys. A: math. Gen 37 (2004)823-839. The authors first transfer the Richards’ partial differential equation into an ordinary differential equation and then they numerically solve the ordinary differential equation. Below is their ordinary DE from pages 835 of their paper which we have solve and graphed using Mathematica. For more information on this readers are referred to the above journal.

 

solution = NDSolve[{3/4 (y[x])^4 y'' [x] + 3 (y[x])^3 (y' [x])^2 - 6/250 (y [x])^7 y' [x] +

          3/10 x y' [x] + 1/10 y [x] Š 0, y [0] == 0.5, y' [0] == 0.1}, y, {x, 0, 15}]

 

Plot[y[x] /. solution, {x, 0, 15}]

 

And the solution is:

 

 

Comments: All the results obtained here were tested by other methods, e.g. integral tables, etc.  and seem to be correct. Various equations such as Error Function, Complementary Error Function and Integral of the Complementary Error Function were obtained from various books such as the book by Beck, J.E., et al., Ozisik, M.L., Carslaw and Jaeger, and some others. These books are mainly cover conduction of heat in solids.