College Math Teaching

October 10, 2011

The Picard Iterates: how they can yield an interval of existence.

Filed under: advanced mathematics, differential equations, numerical methods, numerical solution of differential equations, Taylor Series, uniqueness of solution — collegemathteaching @ 9:33 pm

One of the many good things about my teaching career is that as I teach across the curriculum, I fill in the gaps of my own education.
I got my Ph. D. in topology (low dimensional manifolds; in particular, knot theory) and hadn’t seen much of differential equations beyond my “engineering oriented” undergraduate course.

Therefore, I learned more about existence and uniqueness theorems when I taught differential equations; though I never taught the existence and uniqueness theorems in a class, I learned the proofs just for my own background. In doing so I learned about the Picard iterated integral technique for the first time; how this is used to establish “uniqueness of solution” can be found here.

However I recently discovered (for myself) what thousands of mathematicians already know: the Picard process can be used to yield an interval of existence for a solution for a differential equation, even if we cannot obtain the solution in closed form.

The situation
I assigned my numerical methods class to solve y'= t + y^2 with y(0)=1 and to produce the graph of y(t) from t = 0 to t = 3 .

There is a unique solution to this and the solution is valid so long as the t and y value of the solution curve stays finite; note that \frac{\partial }{\partial y}f(t,y)=2y.

So, is it possible that the y values for this solution become unbounded?

Answer: yes.
What follows are the notes I gave to my class.

Numeric output seems to indicate this, but numeric output is NOT proof.

To find a proof of this, let’s turn to the Picard iteration technique. We
know that the Picard iterates will converge to the unique solution.

y_{0}=1

y_{1}=1+\int_{0}^{t}x+1dx=\frac{1}{2}t^{2}+t+1

y_{2}=1+\int_{0}^{t}x+(\frac{1}{2}x^{2}+x+1)^{2}dx=

y_{2}=\frac{1}{20}t^{5}+\frac{1}{4}t^{4}+\frac{2}{3}t^{3}+\frac{3}{2}t^{2}+t+1

The integrals get pretty ugly around here; I used MATLAB to calculate the
higher order iterates. I’ll show you y_{3}

y_{3}=\frac{49}{60}t^{5}+\frac{13}{12}t^{4}+\frac{4}{3}t^{3}+\frac{3}{2}t^{2}+t+1+O(t^{11})

where O(t^{11}) means assorted polynomial terms from order 6 to 11.

Here is one more:

y_{4}=\frac{17}{12}t^{5}+\frac{17}{12}t^{4}+\frac{4}{3}t^{3}+\frac{3}{2}t^{2}+t+1+O(t^{23})

We notice some patterns developing here. First of all, the coefficient of
the t^{n} term is staying the same for all y_{m} where m\geq n.

That is tedious to prove. But what is easier to show (and sufficient) is
that the coefficients for the t^{n} terms for y_{n} all appear to be
bigger than 1. This is important!

Why? If we can show that this is the case, then our ”limit” solution \sum_{k=0}^{\infty }a_{k}t^{k} will have an interval of convergence less than 1. Why? Substitute t=1 and see that the sum
diverges because the a_{k} not only fail to converge to zero, but they
stay greater than 1.

So, can we prove this general pattern?

YES!

Here is the idea: y_{m}=q(t)+p(t) where p(t) is a polynomial of order m
and q(t) is a polynomial whose terms all have order m+1 or greater.

Now put into the Picard process:

y_{m+1}=1+\int_{0}^{t}((q(x)+p(x))^{2}+xdx=

1+\int_{0}^{t}((q(x)^{2}+2p(x)q(x))dx+\int_{0}^{t}(p(x))^{2}+xdx

Note: all of the terms of y_{m+1} of degree m+1 or higher must come from
the second integral.

Now by induction we can assume that all of the coefficients of the
polynomial p(x) are greater than or equal to one.

When we ”square out” the polynomial, the coefficients of the new
polynomial will consist of the sum of positive numbers, each of which is
greater than 1. For the coefficients of the polynomial (p(x))^{2} of
degree m or higher: if one is interested in the k^{\prime }th
coefficient, one has to add at least k+1 numbers together, each of which
is bigger than one.

Now when one does the integration on these particular terms, one, of course,
divides by k+1 (power rule for integration). But that means that the
coefficient (after integration) is then greater than 1.

Here is a specific example:

Say p(x)=a+bx+cx^{2}+dx^{3}

Now p(x)^{2}=a^{2}+(ab+ab)x+(ac+ca+b^{2})x^{2}+(ad+da+bc+cb)x^{3}+\{O(x^{6})\}

Remember that a,b,c,d are all greater than or equal to one.

Now p(x)^{2}+x=a^{2}+(ab+ab+1)x+(ac+ca+b^{2})x^{2}+(ad+da+bc+cb)x^{3}+\{O(x^{6})\}

Now when we integrate term by term, we get:

\int_{0}^{t}(p(x))^{2}+xdx=a^{2}x+\frac{1}{2}(ab+ab+1)x^{2}+\frac{1}{3}(ac+ca+b^{2})x^{3}+\frac{1}{4}(ad+da+bc+cb)x^{4}+\{O(x^{7})\}

But note that ab+ab+1>2,ac+ca+b^{2}\geq 3, and ad+da+bc+cb\geq 4

Since all of the factors are greater than or equal to 1.

Hence in our new polynomial approximation, the order 4 terms or less all
have coefficients which are greater than or equal to one.

We can make this into a Proposition:

Proposition
Suppose p(x)=\sum_{j=0}^{k}a_{j}x^{j} where each a_{j}\geq 1.

If q(t)=\sum_{j=0}^{2k+1}b_{j}x^{j}=1+\int_{0}^{x}((p(t))^{2}+t)dt

Then for all j\leq k+1,b_{j}\geq 1.

Proof. Of course, b_{0}=1,b_{1}=a_{0}^{2}, and b_{2}=\frac{2a_{0}a_{1}+1}{2}

Let n\leq k+1.

Then we can calculate: (since all of the a_{n-1},a_{n-2},....a_{1} are
defined):

If n is odd, then b_{n}=\frac{1}{n}(2a_{0}a_{n-1}+2a_{1}a_{n-2}+...2a_{\frac{n-3}{2}}a_{\frac{n+1}{2}}+(a_{\frac{n-1}{2}})^{2})\geq \frac{1}{n}(2\ast \frac{n-1}{2}+1)=1

If n is even then b_{n}=\frac{1}{n}(2a_{0}a_{n-1}+2a_{1}a_{n-2}+....2a_{\frac{n-1}{2}}a_{\frac{n+1}{2}})\geq \frac{1}{n}(2\ast \frac{n}{2})=1

The Proposition is proved.

Of course, this possibly fails for b_{n} where n>k+1 as we would fail to
have a sufficient number of terms in our sum.

Now if one wants a challenge, one can modify the above arguments to show that the coefficients of the approximating polynomial never get ”too big”; that is, the coefficient of the k^{\prime }th order term is less than, say, k.

It isn’t hard to show that b_{n}\leq \max a_{i}^{2} where i\in\{0,1,2,...n-1\}

Then one can compare to the derivative of the geometric series to show that
one gets convergence on an interval up to but not including 1.

September 23, 2011

Defeating DFIELD6 (a direction field solver for Matlab)

Filed under: calculus, derivatives, differential equations, mathematics education, uniqueness of solution — collegemathteaching @ 5:01 pm

I was teaching students about the existence and uniqueness theorems for first order differential equations; what I taught was the following version:

Given y'=f(t,y) and an open rectangle R = {(t,y), a<t<b, c < y < d} and an initial point (t_0,y_0) \in R
a) If f(t,y) is continuous over R then there is at least one solution y(t) where y(t_0) = y_0 and the solution is valid so long as the graph of y(t) stays in R .
b) if, in addition, \frac{\partial f}{\partial y} is continuous over R then the solution is unique.

Ok, ok, I know that there are stronger statements of the uniqueness theorem but this is fine for a first pass.

We then talked about what this means: basically, if the conditions of the uniqueness theorem are met, then solution curves for the same differential equation cannot intersect (if they did, this would be a point where uniqueness does not hold).

So, for grins, I plotted the direction field (using Matlab and dfield6) for y' = (ty)^{4/5} and showed how the solver struggled along the y = 0 axis due to the partial derivative not existing there. But then:

Check out the solutions below the y = 0 axis. What is going on there? We should have uniqueness!

Then note:
1. These solutions are NOT following the direction fields…
2. The direction field itself is wrong…do you see why? And
3. Note that the “solution” curves have some local maximums and local minimums…that is impossible!

Let’s look at the differential equation again:
y' = (ty)^{4/5} . We see the following immediately:
1. y' is always greater than or equal to zero and can only equal zero when ty = 0 . Hence the lines with “negative” slope in the the direction field are bogus.
2. Any curve can only have relative maximums or minimums when t = 0 and when y = 0 . Our solver put in some bogus ones.

We can work on this “by hand” and see that the unique solution meeting y(0) = y_0 is y(t) = (\frac{1}{9}t^{\frac{9}{5}} + y_0^{\frac{1}{5}})^5 which is valid so long as y \neq 0

We can also examine the differential equation itself to see that the “curves of constant slope” are of the form ty = \pm c ; example: along the hyperbolas y = 1/t and y = -1/t we should have slopes equal to 1 in the direction field. Clearly this is not being accurately represented by dfield6 for this differential equation!

So what is going on?

Hint: try to solve this differential equations with Matlab Runge-Kutta methods (and others); what you’ll find is that the software is using complex roots for the “raised to the 4-5’ths” function.

Upshot: we are still smarter than the software.

May 16, 2011

Why is 0! = 1?

Filed under: advanced mathematics, calculus, integrals, mathematics education, media, uniqueness of solution — collegemathteaching @ 2:39 pm

This is in response to the following article which says, in part:

I do want to bring up one interesting case study I came across that points in favor of the “math is invented” side of the debate. My friends over at the popular blog Ask a Mathematician, Ask a Physicist did a great post a while ago addressing one of their readers’ questions: What is 0^0?

The reason this question is a head-scratcher is that our rules about how exponents work seem to yield two contradictory answers. On the one hand, we have a rule that zero raised to any power equals zero. But on the other hand, we have a rule that anything raised to the power of zero equals one. So which is it? Does 0^0 = 0 or does 0^0 = 1? […]

Indeed, the Mathematician at AAMAAP confirms, mathematicians in practice act as if 0^0 = 1. But why? Because it’s more convenient, basically. If we let 0^0=0, there are certain important theorems, like the Binomial Theorem, that would need to be rewritten in more complicated and clunky ways.

I provided an answer to the 0^0 question in another post. But the binomial theorem (which indeed starts with (a + b)^0 uses the binomial coefficients; and for each n we eventually have to do with n choose k notation {n \choose k} which involves a 0! in the denominator when k = n or when k = 0 . And yes, the formulas don’t work unless 0! = 1 .

But is that the only reason that 0! = 1 ?

Hardly. As before, we’ll use the properties of analysis to see why this should be true.
Note: I can recommend this internet article for those who are unfamiliar with some of the topics that we are going to discuss; in fact, that article is better than what I am about to write. 🙂

So, can the factorial function n! = n(n-1)(n-2)....(3)(2)(1) be extended so as to be continuous on at least the non-negative part of the real line? The answer is “yes”; the usual extension is the so called gamma function \Gamma (x) = \int^{\infty}_{0} t^{x-1} e^{-t} dt

It is an easy exercise in improper integration, integration by parts and mathematical induction to see that \Gamma (n + 1) = n! for all n \in {1, 2, 3, ....} and that \Gamma (1) = 1 ; hence it makes sense to set 0! =1 .

It is a more advanced exercise to see that the gamma has continuous derivatives of all orders defined for all positive reals (hint: differentiate under the integral sign) and that, in fact, the gamma function is real analytic (has a Taylor series defined on some open interval at every point ).

Strictly speaking, the gamma function is NOT the only possible way to extend the factorial function to an analytic function (example: \Gamma (x) + sin(\pi x) = \Gamma (x) on the non-negative integers ). But if we add the following restrictions:
1. f(1) = 1
2. f(x+1) = xf(x) (a property the factorial function has) and
3. f is logarithmically convex (e. g., ln(f(x) is convex in x then
the Bohr-Mollerup theorem shows that the gamma function is the only extension.

Now of course, the key here is the recursive property (2); this forces the 0! =1 definition.

March 6, 2010

Why We Shouldn’t Take Uniqueness Theorems for Granted (Differential Equations)

Filed under: differential equations, partial differential equations, uniqueness of solution — collegemathteaching @ 11:07 pm

I made up this sheet for my students who are studying partial differential equations for the first time:

Remember all of those ”existence and uniqueness theorems” from ordinary differential equations; that is theorems like: “Given

y^{\prime }=f(t,y) where f is continuous on some rectangle
R=\{a<t<b,c<y<d\} and (t_{0},y_{0})\in R, then we are guaranteed at least one solution where y(t_{0})=y_{0}. Furthermore, if \frac{\partial f}{\partial y} is continuous in R then the solution is unique”.

Or, you learned that solutions to
y^{\prime \prime }+p(t)y^{\prime}+q(t)y=f(t), y(t_{0})=y_{0}, \ y^{\prime}(t_{0})=y_{1} existed and were unique so long as p,q, and f were continuous at t_{0}.

Well, things are very different in the world of partial differential
equations.

We learned that u(x,y)=x^{2}+xy+y^{2} is a solution to xu_{x}+yu_{y}=2u
(this is an easy exercise)

But, can attempt a solution of the form u(x,y)=f(x)g(y).
This separation of variables technique actually works; it is an exercise to see that u(x,y)=x^{r}y^{2-r} is also a solution for all real r!!!

Note that if we wanted to meet some sort of initial condition, say, u(1,1)=3, then u(x,y)=x^{2}+xy+y^{2}, and u(x,y)=3x^{r}y^{2-r} provide an infinite number of solutions to this problem. Note that this is a simple, linear partial differential equation!

Hence, to make any headway at all, we need to restrict ourselves to studying very specific partial differential equations in situations for which we do have some uniqueness theorems.

Implicit Differentiation and Differential Equations

Filed under: calculus, differential equations, implicit differentiation, Uncategorized, uniqueness of solution — Tags: , — collegemathteaching @ 2:08 pm

In my “business calculus class”, we were studying implicit differentiation.
We had a problem:
Find dy/dx if x/y - x^2 = 1
I showed them three ways to do the problem, all of which yield different looking answers:

x - x^2y = y Differentiate both sides: 1 - 2xy -(x^2)dy/dx = dy/dx which yields:
dy/dx = (1-2xy)/(1 + x^2)

Method 2: do directly:
(y - (x)dy/dx)/y^2 -2x= 0  This leads to dy/dx = (y - 2xy^2)/x

Of course that looks different; but we can always solve for y and do it directly:
x/y = 1+ x^2 which yields y = x/(x^2+1) which yields the easy solution:
dy/dx = (1 - x^2)/((x^2 + 1)^2

Now one can check that all three solutions are in fact equal on the domain y \neq 0

But here is the question that came to mind:  in the first two methods we had two different two variable equations: (y - 2xy^2)/x^2 = dy/dx = (1-2xy)/(1 + x^2)

So what does this mean for y ?  Is it uniquely determined?

Answer:  of course it is: what we really have is f(x,y)=g(x,y) = dy/dx whose solution IS uniquely determined on an open rectangle so long as f and g are continuous and  \partial f/y and \partial g/y are continuous also. 🙂

But I didn’t realize that until I took my morning swim. 🙂

This is the value of talking to a friend who knows what he/she is doing: I was reminded that f(x,y)=g(x,y) = dy/dx means that f = dy/dx and g = dy/dx indeed have unique solutions that have the same slope at a common point, but with just this there is no reason that the solutions coincide over a whole interval (at least without some other condition).

So we have something to think about and to explore; I don’t like being wrong but I love having stuff to think about!

Now, of course, we have “different” differential equations with the same solution; yes, there is a theory for that.  I’ve got some reading to do!

Blog at WordPress.com.