Symbolic Calculations with Mathematica.

This demonstration shows some examples of deriving formulas from the class using Mathematica's symbolic functionality. It's important to realize that Mathematica can make a lot of the complicated-looking algebraic and Taylor series manipulations that we do in this subject a lot easier.

Example Problem 1. Derive an O(h^2) error formula for f''(x) using Mathematica.

We start by writing down the Taylor series for f around x. What's the Mathematica command for Taylor series again?

?Series

Ok, we're good. Now we go ahead and work out the Taylor series expansion for f(x+h) and f(x-h).

Clear[p];

p =Series[f[x +h],{h,0,3}]
m  = Series[f[x-h],{h,0,3}]

We now want to combine p and m to leave the f''[x] term. Let's try it:

p + m

We now want to solve this expression for f''[x]. It's not so easy in this form, because Mathematica keeps series expansions in a funny internal form. We can get around this by converting these to regular polynomials with "Normal".

ppm = Normal[p+m]

Now we're ready to make this an equation. In Mathematica, equations are given by the == sign. We try it:

Solve[f[x+h] + f[x-h] == ppm,f''[x]]

If we were really good at Mathematica, we'd probably write all that out in one line:

Solve[f[x+h] + f[x-h] == Normal[Series[f[x+h],{h,0,3}]] + Normal[Series[f[x-h],{h,0,3}]],f''[x]]

Of course, now we can play around with the formula to see if we can get something better.

Solve[f[x+h] + f[x-h] == Normal[Series[f[x+h],{h,0,5}]] + Normal[Series[f[x-h],{h,0,5}]],f''[x]]

As you can see, it's not so much fun to just keep going in the power series, since we just get this f''''(x) term. We can't evaluate _that_! But we can play around with adding some more points. How about f[x + 2h] and f[x-2h]?

Clear[h];
Solve[f[x+2h] + f[x-2h] == Normal[Series[f[x+2h],{h,0,5}]] +
Normal[Series[f[x - 2h], {h,0,5}]],f''[x]]

Score! We have an h^4 f''''(x) term in both expressions. We want to combine them to eliminate this term. To do that, we need to be a little better at eliminating variables from systems of equations with Mathematica. Observe that we really have two equations here:

eq1 := f[x+h] + f[x-h] == Normal[Series[f[x+h],{h,0,5}]] + Normal[Series[f[x-h],{h,0,5}]]

eq2 := f[x+2h] + f[x-2h] == Normal[Series[f[x+2h],{h,0,5}]] +
Normal[Series[f[x - 2h], {h,0,5}]]
eq1 
eq2

We want to eliminate f''''[x] and solve the remaining single equation for f''[x]. That's no problem!

Solve[Eliminate[{eq1,eq2},f''''[x]],f''[x]]

Now this involved a certain amount of algebraic smartitude on our part because we saw immediately how to combine the f[x+h] and f[x-h] and f[x+2h] and f[x-2h] terms to obtain an expression only involving f''''[x] and stuff we wanted. 
But that was really a matter of eliminating f'[x] and f'''[x], wasn't it? We can do _that_ with Mathematica, too!

So here's a way to obtain derivative estimators from Mathematica without any intelligence on our part whatsoever. First, we come up with expansions for f at the gridpoints we'd like to evaluate at:

eqA := f[x+h] == Normal[Series[f[x+h],{h,0,5}]]
eqB := f[x-h] == Normal[Series[f[x-h],{h,0,5}]]
eqC := f[x+2h] == Normal[Series[f[x+2h],{h,0,5}]]
eqD := f[x-2h] == Normal[Series[f[x-2h],{h,0,5}]]
Eqs = {eqA, eqB, eqC, eqD}

Now we start eliminating all of the terms we don't like. In principle, we should only be able to eliminate three variables from a system of four equations. But we get lucky in the algebra because this is so structured:

Eliminate[Eqs,{f'[x],f'''[x],f''''[x],f'''''[x]}]
Solve[%,f''[x]]

What's that %? In Mathematica, the % stands for "the result of the last calculation". Commonly, you read it as "that". So the above reads, "get one equation from our system of four equations by eliminating f', f''', f'''', and f'''''. Now solve THAT for f''."

What's the error in this estimate?

Remember, we took expansions to the 5th order, so the error term must have been at least O[h^6]. If it was actually O[h^7], we could have repeated this process with one more term in the expansion. Can we?

eqA := f[x+h] == Normal[Series[f[x+h],{h,0,6}]];
eqB := f[x-h] == Normal[Series[f[x-h],{h,0,6}]];
eqC := f[x+2h] == Normal[Series[f[x+2h],{h,0,6}]];
eqD := f[x-2h] == Normal[Series[f[x-2h],{h,0,6}]];
Eqs = {eqA, eqB, eqC, eqD};
Eliminate[Eqs,{f'[x],f'''[x],f''''[x],f'''''[x]}];
Solve[%,f''[x]]

Rats. This estimator includes an f'^(6)[x] term which wasn't eliminated.  So the error really _was_ O[h^6]. 

Still, this Mathematica stuff is really cool, isn't it?