Mathematica Lesson II

Physics 80301 Fall 2005

Solving Algebraic Equations

Single Nonlinear Equation

Let us consider a nonlinear equation such as (Sqrt[4 - (x^2)])= tan(x) that occurs in the quantum mechanical square-well problem. As a first step in finding the solutions, plot the left and right sides of the equality and look for intersections. In this case, as seen below, there is an intersection near x = 1.06. (To find intersection hold down control key and move cursor over intersection point of the graph.)

Plot[{(Sqrt[4 - (x^2)]), Tan[x]}, {x, 0, 2}, 
      PlotStyle -> {{Thickness[0.01], 
            RGBColor[1, 0, 0]}, {Thickness[0.01], RGBColor[0, 0, 0]}}, 
      PlotRange -> {0, 2}, AspectRatio -> Automatic]

To get an accurate value use FindRoot[(Sqrt[4 - (x^2)]) - Tan[x],{x,1.04}], where the value 1.04 is the approximation from the plot.
This routine uses Newton's method starting from the approximate root and homes in on the exact root.

FindRoot[(Sqrt[4 - (x^2)]) - Tan[x], {x, 1.04}]

Let's check the root. As we see below, the equation is solved to parts in (10^16).

(Sqrt[4 - (x^2)]) - Tan[x] /.  %

Several Polynomial Equations

As an example consider the pair of equations (x^2)+(y^2)=4 and x+y=1. These equations can be solved as follows:

Solve[{(x^2)+(y^2)==4, x + y == 1},{x, y}]

Note the use of == above to designate equality.
To check the solution, we substitute back into the equations

Simplify[{(x^2)+(y^2),x + y}/. %]

To get a numerical solution, use NSolve[]  instead of Solve[]:

NSolve[{(x^2)+(y^2)==4, x + y == 1},{x, y}]

Check the solutions:

Simplify[{(x^2)+(y^2),x + y}/. %]

Several Linear Equations

Use matrix algebra!
 A . x = b  is solved as x = Inverse[A] . b

A  =  {{1, 3, 2}, {-3, 1, 2}, {-2, -1, 3}};  
    MatrixForm[A]

b  =  {11/2, 1/2, 2}

x  =  Inverse[A ] .  b

Check the solution

A  .  x  - b

Solving Differential Equations

Analytical Solutions

Consider as an example the damped oscillator equation x''(t) +(3/2) x'(t) + 4 x(t) = 0.

Clear[x];
sol  = 
      DSolve[{x''[t] +  (3/2) x'[t] +  4  x[t] == 0, 
          x[0] == 1, x'[0] == 0}, x[t], t]

Plot[x[t] /. sol, {t, 0, 6}, PlotRange -> {- .5, 1}, 
      PlotStyle -> {RGBColor[1, 0, 0], Thickness[0.01]}, 
      AxesLabel -> {time , Amplitude}]

Numerical Solutions

Consider the equation for a finite-amplitude pendulum: θ''(t) + (ω^2)sin(θ) = 0. Initial conditions:  θ(0) = (Subscript[θ, 0]), and θ'(0) = 0.  (Set ω=1)

NDSolve[{θ''[t] +  Sin[θ[t]] == 
          0, θ[0] == 15  Pi/16, θ'[0] == 
          0}, {θ[t]}, {t, 0, 7.4}]

Plot[θ[t] /. %, {t, 0, 7.4}, 
      PlotStyle -> {RGBColor[1, 0, 1], Thickness[0.015]}, 
      AxesLabel -> {time , Angle}]

Now let's consider an alternative way of plotting the results. As a first step define the function θ[t].

u  =  θ[t]  /.   First[%%]

tbl  =  
        N[Table[{Sin[u], -Cos[u]}, {t, 0, 7.4, 0.2}], 3];

ListPlot[tbl, AspectRatio -> Automatic, 
      PlotStyle -> {RGBColor[1, 0, 0], PointSize[0.03]}, 
      AxesLabel -> {X, Y}]

Least Squares Approximation

As an example, we fit a straight line y = a(x+2) to the following data points

points = {{0, -1}, {1, 2}, {3, 7}, {4, 6}, {7, 9}, {11, 
            14}};

We use the option  DisplayFunction->Identity to suppress temporarily the plot.

plotpoints  = 
        ListPlot[points, PlotStyle -> PointSize[0.03], 
          DisplayFunction -> Identity];

We fit the points to x+2.  For a quadratic fit one would replace x+2 by the array  {1,x,(x^2)} .

f  =  Fit[points, x + 2, x]

Now we set up a plot of the best fit line and suppress it.

plotline = 
        Plot[f, {x, -3, 14}, 
          PlotStyle -> {RGBColor[1, 0, 0], Thickness[0.01]}, 
          DisplayFunction -> Identity];

Plot the points and line together using Show[].  We restore DisplayFunction to its default value in the plots
by using the option DisplayFunction->$DisplayFunction

Show[plotpoints, plotline, AxesLabel -> {X, Y}, 
      DisplayFunction -> $DisplayFunction]

Power Series

As an example consider f(x) = (e^(sin(x))).  The following gives a series expansion about point x=0  with remainder
of  O(([x])^8).

Series[(E^Sin)[x],{x, 0, 7}]

To convert the series to a polynomial use Normal[]

poly  =  Normal[%]

Compare in the range [-1,1]

Plot[{(E^Sin)[x], poly}, {x, -2, 2}, 
      PlotStyle -> {{RGBColor[0, 0, 0], 
            Thickness[0.01]}, {RGBColor[1, 0, 0], 
            Thickness[0.01]}}]