Differential Equations I



Integration of Ordinary Differential Equations

Press - Chapter 15 first edition, Chapter 16 second edition

Note the section headings -- Runge-Kutta, Adaptive-Step size Control, Modified Midpoint Method, Richardson Extrapolation and ..., Second-Order Conservative Equations, Stiff sets of Equations.

Euler’s method is so simple that you need only remember the name.

Runge-Kutta is a method of using a succession of points to estimate derivatives. It enables one to get convergence of order h5-6 without actually taking any derivatives. It works for mathematical functions that are easy to calculate. It should be used for calculating very accurate estimates of the types of functions that are found to 10 digits in tables. The adaptive size control increases the efficiency of an already efficient method but the changing of step sizes introduces errors and problems.

The modified midpoint method is the closest Press comes to my favorite scheme (sketched at the end of this lecture), which is mid-point with predictor corrector. Richardson extrapolation can be used to increase the nominal convergence rate from h2 to h4 ..\integration\Richardson.htm, My recommendation is to wait until the end rather than trying to do it intermediately as in Press. Second-order conservative equations are our stock in trade and the Numerov method is very frequently used to reduce the number of operations when the first derivative is not wanted. I usually prefer the extra information given by also calculating the first derivative separately.

If you are attempting to solve the equation numerically, it is probably stiff. Press’s section on Stiff differential equations is must reading. The predictor corrector part of mid point trap is the beginnings of a solution to this problem. The Schroedinger equation in one dimension is a stiff differential equation. If you start at r=0 in the hydrogen atom ground state equation and integrate outward with 16 digit precision it blows up by r=16 aB. This will be discussed in much more detail.

Following Press sort of (integration of Ordinary Differential Equations)

The standard equations of physics are second order differential equations of the form

[pic] Eqn 1

There are three functions of x here

These can be rewritten as

[pic]Eqn 2

The function f (x,y(x),z(x)) is explicitly given in terms of y(x) and z(x). The function z(x) is the integral

[pic]

The function y(x) is

[pic]

This reduces a second order differential equation to two simultaneous first order differential equations that are to be solved together. Note that if there were no y or z in f that the solution for z(t) would be[pic]. The reason that this is a differential equation rather than simple quadrature is due to that fact that it is necessary to find y and z before evaluating f.

Deviating from Press Euler’s method

The simplest method for solving these equations is Euler’s method.

x=x0

y=y0

z=z0

h=.001 (think small)

5 y=y+h*z error of order h2

z=z0+h*(r(x,y0,z)-q(x,y0,z)) error of order h2

y0=y

write(1,*)x,y,z

x=x+h

if(x.lt.x1)goto 5

stop

At each step, an error of h2 is made so that after 1/h steps, the total error is of order h. One can plot this and extrapolate it out

[pic] Eqn 3

Solve for y(0).

[pic]Eqn 4

For the case shown h1=2 h2 so that

[pic] Eqn 5

This is to say that the correct answer moves away from the second answer by the difference between the two answers which is shown as e in the above figure.

The values at the beginning of the interval are the only ones used to estimate the entire interval. In terms of integrals we have said.

[pic] Eqn 6

[pic] Eqn 7

Note that I have made no attempt to use the final y in finding z in the above code. The method above which is the simplest possible, is the recommended method for the first pass. If it extrapolates to an answer that you want to publish you will need to remember that it is called Euler’s method and be sure to mention that you extrapolated to arrive at the final answer, otherwise the referee will understand how crudely you solved the problem and object.

A bit more analytic accuracy

The integrand f can be expanded to one higher order as

[pic] Eqn 8

and likewise the integrand z(x’) can be expanded

[pic] Eqn 9

this in turn can be integrated to yield

[pic] Eqn 10

Then

[pic]

[pic] Eqn 11

This actually works very well, if someone is willing to take the partials of the function f with respect to its variables. ..\WaveFunction\1-d\DiffEq\SeqnIn1d.docx. Note that in general evaluating the partials on the computer takes no more time than evaluating f itself making the method above two to three times faster than the Runge-Kutta algorithm that we are about to derive.

Runge-Kutta

The notion is to evaluate at different values of x, y, and z in the interval between x0 and x, and to make a weighted average of these evaluations expand to be the same as the analytic expansion given above. This allows us to get the benefit of the derivative expansion without actually taking derivatives.

[pic]Eqn 12

which will equal the analytic expansion for a=b=h/2; p=q=s=1.

It will also work for a=0,b=1,p=q=s=1/2.

[pic] Eqn 13

In the case a=b=h/2, p=q=s=1, we have (=(=h/2.

In the case a=0,b=h,p=q=s=1/2, we have (=0,(=h. The first case amounts to end point sampling while the second is midpoint sampling.

By using more points and more algebra, the method can be extended to higher and higher powers of h. The method coded in Press evaluates f at the end points and a group of mid points to get an error of order h5. Unfortunately he does only a single equation.

Stiff differential equations

Consider the rather simple equation

[pic] Eqn 14

with c>0. The analytic solution to this is

[pic] Eqn 15

as can easily be seen by taking the derivative. Euler’s method yields

[pic] Eqn 16

Here we see the origin of an oscillation. If a part of the problem is such that (1-ch) is less than -1, the continual multiplying by this will produce an oscillating solution as part of the difference equations. Now to show how very slight changes can have very large results, use the advanced value of y as the derivative, rather than the current value so that.

[pic] Eqn 17

this becomes

[pic] Eqn 18

or

[pic] Eqn 19

which for c>0 is stable for all step sizes. Note that c ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download