Introduction to Mathematica's
Linear Algebra, Calculus, and Differential Equations Tools
W. Craig Carter MIT IAP 2006
http://pruffle.mit.edu/

Linear Algebra:
There are two "flavors" of applications of linear algebra in Mathematica.  The most simple is the generalization of Solve to several linear equations and variables.  Some related functions are Reduce, Resolve, Eliminate, and FindInstance.

First, one must become familiar with Lists. Lists are collections of other things that can be identified with a symbol:

SantasList  =  List[Naughty,  Nice]

Curly brackets are a shorthand and can be used instead of typing List

\(\(Naughty = \ 
        List["George W.", \ "Tom D.", \ "J. Abramoff"];\)
\)
Nice  =  {Girls,  
          Boys};

Boys =  {"Jimmy",  "Billy"};

Girls  =  {"Di",  "Nancy",  "Marie Claire"};

SantasList

To extract elements of a list

Part[SantasList, 2]

Part[SantasList, 2, 1, 3]

Square brackets are the shorthand for Part, the following are equivalent to Part[SantasList,2,1,3]

SantasList[[2]][[1]][[3]]

SantasList[[2, 1, 3]]

There are lots of functions designed to operate on lists and list structures, Last, Sort, Reverse, Select, etc...
Let's look at some examples:

Length[SantasList]

AList  =  {a,  b, 2, 7, 9,  1.3,  (π/2),  0}

Cos[AList]

Attributes[Cos]

MyFunction[AList]

SetAttributes[MyFunction, Listable]

MyFunction[AList]

Sort[AList]

Reverse[Sort[Select[AList, NumberQ]]]

Select[AList, EvenQ]

Select[AList, PrimeQ]

Perms  =  Permutations[Select[AList, ExactNumberQ]]

MatrixForm[Perms]

MatrixForm[Transpose[Perms]]

An Array is a  List of Lists of equalsize---it is a "hyper-rectangular list"  (Those of you who know what a tensor is---tensors have this form). Perms is an Array:

Dimensions[Perms]

Your Turn to Try (two problems):

1. Look up Table[] and Prime[] in the help browser.

Make a 100 element list of r[i] where
	
	r[i] = Prime[i + 1]/Prime[i]
Plot it.

2. Construct an array, MyArray, with Dimensions[] = {6,6,6} where
	
		MyArray[[i,j,k]] = Sin[((i Pi)/6)] Sin[((j Pi)/6)] Sin[((k Pi)/6)]
Write a function to generalize this last part to dimensions {N,M,P}

Linear Systems of equations and variables

Consider the set of equations
	 x + 2y +  z  +  t  = a
	-x + 4y - 2z         = b
	 x + 3y + 4z + 5t = c
	 x          +  z  +   t = d
We illustrate several different methods to work with these equations and their solutions.	  The methods are meant to illustrate different ways to do the same thing as a way of exploring Mathematica, and not necessarily the best way to proceed.

Start with the matrix of coefficients of the variables, mymatrix:

mymatrix  =  {
{1, 2, 1, 
            1}, 
{-1, 4, -2, 
            0}, 
{1, 3, 4, 5}, 
{1, 0, 1, 
            1}};
mymatrix // MatrixForm

The system of equations will only have a unique solution if the determinant of mymatrix is nonzero.

Det[mymatrix]

Now define vectors for (Overscript[x, ⇀])  and (Overscript[b, ⇀]) in \(A\+_\ x\&⇀\  = \ b\&⇀\)

myx  = {x, y, z, t};

myb  =  {a, b, c, d};

The left-hand side of the first equation will be

(mymatrix . myx)[[1]]

and the left-hand side of all four equations will be

lhs  =  mymatrix . myx;
lhs // MatrixForm

Now define an indexed variable linsys with four entries, each being one of the equations in the system of interest:

linsys[i_Integer]  :=  lhs[[i]]  ==  myb[[i]]

linsys[2]

So the list of equations is:

equationlist  =  Table[linsys[i], {i, Length[myb]}]

Solving the set of equations for the unknowns (Overscript[x,  ->])

linsol  = Solve[equationlist, myx]

Notice that the result is a List of List of rules, and not necessarily solutions---to get the solutions, one has to apply replacement

The next line applies all the lists in the list of rules, so the result is a list of a list:

myx /. linsol

If, one wanted only the first rule (in this linear case, the only solutions ruleset)

myx /. linsol[[1]]

And yet another way, based on \(TraditionalForm\`x\&⇀ = \(\(A\+_\^\(-1\)\) 
            A\+_\ x\&⇀\  = \ \(A\+_\^\(-1\)\) 
            b\&⇀\)\)

inversemymatrix =  Inverse[mymatrix];
inversemymatrix // MatrixForm

inversemymatrix . myb // MatrixForm

Doing the same thing a different way, using Mathematica's LinearSolve function:

?LinearSolve

LinearSolve[mymatrix, myb]

Other matrix operations

When determinants are zero

Create a matrix with one row as a linear combination of the others

myzeromatrix = 
{mymatrix[[1]], 
mymatrix[[2]], 
p*
              mymatrix[[1]]  +  q*mymatrix[[2]]  +  
            r*mymatrix[[4]], 
mymatrix[[4]]};
myzeromatrix // MatrixForm

Det[myzeromatrix]

LinearSolve[myzeromatrix, myb]

This was not expected to have a solution because one of four equations in the system was a linear combination of others in the system. Effectively, we were asking Mathematica to solve a system of three equations in four unknowns. The rank of a square matrix of coefficients is equal to the number of linearly independent equations in the system. The null space of the matrix will be empty when the equations are all linearly independent.

MatrixRank[mymatrix]
MatrixRank[myzeromatrix]

NullSpace[mymatrix]
NullSpace[myzeromatrix]

Try solving this inhomogeneous system of equations using Solve:

zerolhs =  myzeromatrix . myx

zerolinsys[i_Integer]  :=  
      zerolhs[[i]]  ==  myb[[i]]

Table[zerolinsys[i], {i, 4}] // MatrixForm

zerolinsol  = Solve[Table[zerolinsys[i], {i, 4}], myx]

No solution, as expected,  Let's see what happens if we ask Mathematica to solve the homogeneous problem:

zerolinsol  = 
      Solve[Table[
          zerolinsys[i] /. {a -> 0, b -> 0, c -> 0, 
              d -> 0}, {i, 4}], myx]

In this case, Mathematica gives a relationship between the variables, but because there are fewer equations than variables, there is still no unique solution.

To obtain the reduced linearly independent system, one could use  the NullSpace or here we  use Eliminate

Eliminate[Table[zerolinsys[i], {i, 4}], {x}]

Your Turn to Try:

For a state of plane stress, the two-dimensional linear elastic stresses and strains are related by

(({{(Subscript[ϵ, xx])}, {({{(Subscript[ϵ, yy])}, {(Subscript[ϵ, xy])}})}}))=(1/E) (({{1, -ν}, {({{-ν}, {0}}), ({{1}, {0}})}})  ({{({{0}, {0}})}, {2 (1 + ν)}}))(({{(Subscript[σ, xx])}, {({{(Subscript[σ, yy])}, {(Subscript[σ, xy])}})}}))

E is called Young's modulus. It is large for hard materials like diamond and small for something like a slinky.  Because E is interpreted as e in Mathematica, you might want to use a different symbol.  ν is Poisson's ratio and it tells you how much something would "squeeze" out as you step on it: rubber has a Poisson's ratio that is about 1/2 and steel is closer to 1/3.

1. Solve for the stresses ((Subscript[σ, xx]), etc) in terms of the strain.

2. Find the conditions on stress, such that all of the strains are equal

3. Look up Eigenvalues, what are the eigenvalues of the above matrix?

Calculus

Derivatives, Maxima, and Minima

The derivative operator is D[]

Clear[f];
D[f[x], x]

D[f[x], {x, 2}]

D[f[x], y]

D[f[x, y, z], {x, 2}, {y, 3}, z]

Clear[x, y, z];
{x, y, z}  =  {a[t], b[t],  c[t]};
D[f[x, y, z], {x, 2}, {y, 3}, z]

D[f[x, y, z], t]

Let's do an example:

f[x_, 1]  :=  x;
f[x_,  n_]  :=  (((2 f[x, n - 1]))^f)[x, n - 1]

OurLuckyDay  =  f[t, 4]

OurLuckyDay /. t -> 0

Limit[OurLuckyDay, t ->  0]

Limit[D[OurLuckyDay, t], t -> 0]

D[OurLuckyDay, t] /. t -> 1 // N

Plot[OurLuckyDay, {t, 0, 1}]

Plot[OurLuckyDay, {t, 0, 1}, PlotRange -> {0, 16}]

dft =  Simplify[D[OurLuckyDay, t]]

Solve[dft == 0, t]

troot  =  FindRoot[dft, {t,  .1,  .5}]

OurLuckyDay /. troot

FindMinimum[OurLuckyDay, {t, 1}]

FindMinimum[OurLuckyDay, {t,  .1,  .5}]

NMinimize[{OurLuckyDay, t >  0}, t]

Integration

Clear["Global`*"];

AMessyExpression  =  Log[x Sin[x]]/((1/x))

DMess  =  D[AMessyExpression, x]

Integrate[DMess, x]

Notice that this next indefinite integral takes a while longer than the indefinite integral

Integrate[DMess, {x, 0, E}] // N

(AMessyExpression /. x -> e)  -  (AMessyExpression /. 
          x -> 0)

(AMessyExpression /. x -> E)  -  
      Limit[AMessyExpression, x -> 0]

N[(AMessyExpression /. x -> E)  -  
        Limit[AMessyExpression, x -> 0]]

Some integrals don't exist in closed form:

Integrate[Sin[x]/Sqrt[(x^2  +  a^2)], x]

But, definite integrals exist in certain cases.

UglyInfiniteIntegral  =  
      Integrate[
        Sin[x]/Sqrt[(x^2  +  a^2)], {x, 0,  ∞}]

Assumptions should point to a list, either connected by logical ands (i.e., And[a,b] is the same as a&&b which is true only if both a and b are true).

UglyInfiniteIntegral  =  
      Integrate[Sin[x]/Sqrt[(x^2  +  a^2)], {x, 0,  ∞}, 
        Assumptions -> Re[a^2] > 0]

Numerical integration is also possible

NIntegrate[Sin[x]/Sqrt[(x^2  +  2^2)], {x, 0, 10}]

Your Turn to Try:

Write a function Snint(z) that calculates:
		
		Snint(z) = (∫_0^z) sin[
            2 π x sin (2  πx)] dx
Plot Snint(z) with respect to z from z=0 to z=3.

Ordinary  Differential Equations

ODE's combine linear algebra and calculus, one asks to solve sets of equations that contain derivatives, for example

MyODE  =  
      D[y[t], {t, 2}]  -  t  D[y[t], t]  +  y[t]  ==  0

DSolve[MyODE, y[t], t]

This has two constants, called C[1] and C[2], that indicate that two boundary conditions are needed for a particular solution

MyParticularSolution  =  
      DSolve[{MyODE, y[0] == 1,  y'[0] == 0}, y[t], t]

y[t] /. MyParticularSolution[[1]]

Plot[y[t] /. MyParticularSolution[[1]], {t, 0, 2}]

Some ODE's do not have closed form solutions

AnotherODE  =  
      D[y[t], {t, 2}]  - Sin[t]  D[y[t], t]  +  y[t]  ==  0

DSolve[{AnotherODE, y[0] == 1,  y'[0] == 0}, y[t], 
      t]

In which case, a numerical version of DSolve, called NDSolve, can be used to find an "Interpolation Function" (which is a recipe to find solutions at various values)

NumericalInterpolator =  
      NDSolve[{AnotherODE, y[0] == 1,  y'[0] == 0}, 
        y[t], {t, 0, 50}]

Plot[y[t] /. NumericalInterpolator[[1]], {t, 0, 50}]

Plot[y[t] /. NumericalInterpolator[[1]], {t, 0, 55}]

Your Turn to Try:

Consider steady-state diffusion (say of heat or a chemical species) across a finite block with diffusivity D.  This is a one-dimensional problem:

	(d/(d x))(D ((d T)/dx))  =  
        0

       
and suppose the boundary conditons are T(x=0) = 1, T(x=1) = 0

Suppose the diffusivity of the block is given by D =  1 + (e^(-10 (x^2))).

1. Calculate a closed form approximation
to T(x) by replacing D with its spatial average <D> = (∫_0^1) D (x) dx

2. Calculate a numerical approximation to T(x).