:

## Two Illustrations of "Is g the same as f?"

Here are two illustrations for how one might want to check to see if g is the same as f. The attached file is a Maple 10 worksheet.

=== Content of Worksheet follows [MaplePrimes moderator]

For the original question, see this forum posting.

Hello, Tom. Newbie? I would hate to pass myself off as an old hand in the august company of these Maple experts! On the other hand, I have been doing this a long time, haven't I?

Here is what my question was all about. I will give an easy example first.

Suppose a newbie (even to me) is asked to solve the simple ordinary differential equation y' + y = 0 with y(0) = 1. Being a real newbie, he might say the following: Ha! I know the solution for this equation. I will define it here.

> y:=t->exp(-3.14*t);

But, having been admonished to always check his work, he does this:

> diff(y(t),t)+Pi*y(t);

Of course, a newbie that is an aspiring engineer might say that this is industrial grade zero.

Upon seeing this response, among other things, I would ask that the graph of this last be drawn.

> plot(-3.14*exp(-3.14*t)+Pi*exp(-3.14*t),t=0..1);

At least now, there is an estimate of how much error there is over the interval [0,1].

> restart;

A year or so later, and with much more experience, this potential engineer takes my PDE class. First thing you know, there is a discussion of the heat equation on a disk. If the solutions are independent of , then the partial differential equation to be solved reduces to = .

Bessel's equation comes up. The student might go looking for the zeros of the Bessel functions. (From the perspective of the mathematics involved, this is comparable to approximating with 3.14 in the above simple problem.) For simplicity, let's compute only 10 of the Bessel zeros.

> for n from 1 to 10 do
zero[n]:=evalf(BesselJZeros(0,n)):
od:
n:='n':

If the initial distribution was , then Fourier coefficients would be computed.

> for n from 1 to 10 do
A[n]:=int((1-r^2)*BesselJ(0,zero[n]*r)*r,r=0..1)/
int(BesselJ(0,zero[n]*r)^2*r,r=0..1):
od:
n:='n':

Then, a potential solution for the equation would be created.

> u:=(t,r)->sum(A[n]*BesselJ(0,zero[n]*r)*exp(-zero[n]^2*t),n=1..10):

Having been told to check and check and check solutions, the student does.
First, the initial value can be checked. HERE IS THE ISSUE: is the same as u(0,r)?

I draw graphs, off-setting u(0,r) a little so that you can tell them apart.

> plot([1-r^2,u(0,r)+0.01],r=0..1,color=[red,green]);

HERE IS THE ISSUE, again. Is the following difference the same as zero. This is a check of the PDE.

> simplify(diff(r*diff(u(t,r),r),r)/r-diff(u(t,r),t));

I think no amount of simply, combine, collect, ... will make that zero and I think I know why. So, I draw a graph of this difference.

> plot3d(%,r=0..1,t=0..1,axes=normal,orientation=[-35,70]);

Even I agree that this graph is a graph of industrial grade zero!

Of course, it the student had used the exact zeros of the Bessel functions, things would have worked differently.

> v:=(t,r)->sum(A[n]*BesselJ(0,BesselJZeros(0,n)*r)*
exp(-BesselJZeros(0,n)^2*t),n=1..10);

> simplify(diff(r*diff(v(t,r),r),r)/r-diff(v(t,r),t));

And, the graph of the solution could be drawn and animated.

> plots[animate3d]([r,theta,v(t,r)],r=0..1,theta=-Pi..Pi,t=0..1,
coords=cylindrical,axes=normal,orientation=[45,65],
frames=25);

This post generated using the online HTML conversion tool 