y(x + h) - Imperial College London
Download
Report
Transcript y(x + h) - Imperial College London
Numerical Solution of
Differential Equations
3rd year JMC group project ● Summer Term 2004
Supervisor: Prof. Jeff Cash
Saeed Amen
Paul Bilokon
Adam Brinley Codd
Minal Fofaria
Tejas Shah
© Imperial College London
Agenda
Adam: Differential equations and numerical
methods
Paul: Basic methods (FE, BE, TR). Problem:
“circle”
Saeed: Advanced methods (4th order).
Problem: Kepler’s equation
Minal: Stiff equations
Tejas: Derivation of 6th order method
© Imperial College London
Agenda
Adam: Differential equations and numerical
methods
Paul: Basic methods (FE, BE, TR). Problem:
“circle”
Saeed: Advanced methods (4th order).
Problem: Kepler’s equation
Minal: Stiff equations
Tejas: Derivation of 6th order method
© Imperial College London
Introduction to Differential
Equations
Equations involving a function y of x, and its
derivatives
Model real world systems
General Equation:
© Imperial College London
Introduction to Differential
Equations
Simple example
Solution obtained by integrating both sides
Initial values can determine c
© Imperial College London
Introduction to Differential
Equations
Special case when equations do not involve x
E.g.
Initial values y(0) = 1 and y'(0) = 0, solution is
© Imperial College London
Introduction to Differential
Equations
Kepler’s Equations of Planetary Motion
Difficult to solve analytically – we use
numerical methods instead
© Imperial College London
Introduction to Differential
Equations
Forward Euler and Backward Euler
Trapezium Rule
General method
–
–
–
Use initial value of y at x=0
Calculate next value (x=h for small h) using gradient
Call this y1 and repeat
Need formula for yn+1 in terms of yn
© Imperial College London
Agenda
Adam: Differential equations and numerical
methods
Paul: Basic methods (FE, BE, TR). Problem:
“circle”
Saeed: Advanced methods (4th order).
Problem: Kepler’s equation
Minal: Stiff equations
Tejas: Derivation of 6th order method
© Imperial College London
Agenda
Adam: Differential equations and numerical
methods P
Paul: Basic methods (FE, BE, TR). Problem:
“circle”
Saeed: Advanced methods (4th order).
Problem: Kepler’s equation
Minal: Stiff equations
Tejas: Derivation of 6th order method
© Imperial College London
Euler’s Methods
Simplest way of solving
ODEs numerically
Does not always
produce reasonable
solutions
Forward Euler (explicit)
and Backward Euler
(implicit)
© Imperial College London
Forward Euler
y
True solution
(x2, y2)
y2
Slope f(x1, y1)
Slope f(x0, y0)
Approximated
solution
y1
y(x + h) = y(x) + h y’(x)
Explicit
Eloc = ½h2y(2)() h2
First order
Asymmetric
(x1, y1)
y0
(x0, y0)
h
x0
h
x1
x2
x
© Imperial College London
Backward Euler
y
(x2, y2)
y2
True solution
Slope f(x2, y2)
Slope f(x1, y1)
Approximated
solution
(x1, y1)
y1
y0
y(x + h) = y(x) + h y’(x +
h)
Implicit
Eloc = - ½h2y(2)() h2
First order
Asymmetric
(x0, y0)
h
x0
h
x1
x2
x
© Imperial College London
Trapezium Rule
y(x + h) = y(x) + ½ h [y’(x) + y’(x + h)]
Average of FE and BE
Implicit
Eloc = - (h3/12) f(2)() h3
Second order
Symmetric
© Imperial College London
Example Problem: “Circle”
Consider y’’ = -y, with
initial conditions
–
–
y(0) = 1
y’(0) = 0
The analytical solution is
y = cos x, that can be
used for comparison with
numerical solutions
© Imperial College London
Plots for the True Solution
Time series plot (xy)
© Imperial College London
Phase plane plot (zy)
Time Series
y
x
© Imperial College London
Phase Plane Plots: FE & BE
Forward Euler
© Imperial College London
Backward Euler
Phase Plane Plots: TR
FE & BE fail: non-periodic
TR: OK, periodic soln
Why?
© Imperial College London
Symmetricity
TR is symmetric, whereas FE & BE are not
y(x + h) = y(x) + ½ h [y’(x) + y’(x + h)]
h -h
y(x - h) = y(x) - ½ h [y’(x) + y’(x - h)]
X := x - h
y(X + h) = y(X) + ½ h [y’(X) + y’(X +
h)]
© Imperial College London
Symmetricity
TR is symmetric, whereas FE & BE are not
y(x + h) = y(x) + h y’(x)
h -h
y(x - h) = y(x) - h y’(x)
X := x - h
y(X + h) = y(X) + h y’(X + h)
© Imperial College London
Time Step
h = 10 -1
© Imperial College London
h = 10 -2
Agenda
Adam: Differential equations and numerical
methods P
Paul: Basic methods (FE, BE, TR). Problem:
“circle”
Saeed: Advanced methods (4th order).
Problem: Kepler’s equation
Minal: Stiff equations
Tejas: Derivation of 6th order method
© Imperial College London
Agenda
Adam: Differential equations and numerical
methods P
Paul: Basic methods (FE, BE, TR). Problem:
“circle” P
Saeed: Advanced methods (4th order).
Problem: Kepler’s equation
Minal: Stiff equations
Tejas: Derivation of 6th order method
© Imperial College London
After Euler and TR
Can create higher order methods, which have far
smaller global errors
Methods are more complex and require more
computation on each step
But for same step size more accurate
Introduce concept of half step
© Imperial College London
Two Fourth Order Methods
yn+1 – yn = yn + h/2(y’n) + (h2/12)(4y’’n+½+ 2y’’n)
y’n+1 – y’n = h/6(y’’n+1 + 4y’’n+½+ y’’n) (*)
Then to calculate the half-step we use either A or B
–
–
A) yn+ ½ = ½(yn+1 + yn) – h2/48(y’’n+1 + 4y’’ n+½ + y’’n)
B) yn+ ½ = yn + ½y’n – h2/192(-2y’’n+1 + 12y’’ n+½ + 14y’’n)
Which one is more accurate? (next slide)
Solve iteratively and then apply solution to find
derivative (*)
© Imperial College London
Creating Method A
Start with
–
–
–
–
yn+1 – yn = h/6[y’n+1 + 4y’n+½ + y’n] (1)
diff. y’n+1 – y’n = h/6[y’’n+1 + 4y’’n+½ + y’’n] (2)
yn+½= ½(yn+1 + y’n) – h/8(y’n+1 – y’n) (3)
diff. y’n+½= ½(y’n+1 + y’’n) – h/8(y’’n+1 – y’’n) (4)
subs (4) into (1) (to eliminate y’n+½ ) and
eliminate y’n+1 using (2),
then use (2) in (3) to get half-step
© Imperial College London
Comparing Fourth Order Methods
Comparing errors when
solving circle problem
Both A and B produce a
much smaller order of
error than Euler’s
© Imperial College London
Finding Earth’s Orbit Around Sun
Kepler’s Equation
Can use it to find the orbit of planets
Use to find orbit of earth around the sun
Work in two dimensions z and y
z’’ = -(GM z) / (y2 + z2)3/2
y’’ = -(GM y) / (y2 + z2)3/2
Constants and initial conditions in report
© Imperial College London
Results Plot
Solution uses small step
size h = 0.01
Becomes difficult to tell
difference between
methods visually
Using h = 0.1 difference
is more marked
© Imperial College London
Agenda
Adam: Differential equations and numerical
methods P
Paul: Basic methods (FE, BE, TR). Problem:
“circle” P
Saeed: Advanced methods (4th order).
Problem: Kepler’s equation
Minal: Stiff equations
Tejas: Derivation of 6th order method
© Imperial College London
Agenda
Adam: Differential equations and numerical
methods P
Paul: Basic methods (FE, BE, TR). Problem:
“circle” P
Saeed: Advanced methods (4th order).
Problem: Kepler’s equation P
Minal: Stiff equations
Tejas: Derivation of 6th order method
© Imperial College London
Stiff Equations
Certain systems of ODEs are classified as stiff
A system of ODEs is stiff if there are two or
more very different scales of the independent
variable on which the dependent variables are
changing
Some of the methods used to find numerical
solutions fail to obtain the required solution
© Imperial College London
Example
Consider the equation:
2
d y
dx
2
( 1)
dy
y 0
dx
With initial conditions:
y (0) 1
y ' (0) 1
There are two solutions to this problem
© Imperial College London
Example continued
Analytical solution: y=e-x
Unwanted solution: y=e-λx
We will now show how forward Euler is not
stable when solving this problem under certain
circumstances
© Imperial College London
Example continued
Let λ = 103 and let (the step size) h = 0.1
© Imperial College London
Increasing the range of the x-axis
© Imperial College London
When h = 0.001
© Imperial College London
Agenda
Adam: Differential equations and numerical
methods P
Paul: Basic methods (FE, BE, TR). Problem:
“circle” P
Saeed: Advanced methods (4th order).
Problem: Kepler’s equation P
Minal: Stiff equations
Tejas: Derivation of 6th order method
© Imperial College London
Agenda
Adam: Differential equations and numerical
methods P
Paul: Basic methods (FE, BE, TR). Problem:
“circle” P
Saeed: Advanced methods (4th order).
Problem: Kepler’s equation P
Minal: Stiff equations P
Tejas: Derivation of 6th order method
© Imperial College London
A Sixth Order Numerical Method
Has an error term with smallest h degree as h7
Idea is to find values for α, A, B, C, D, E, F, C,
D, E, F such that:
© Imperial College London
Derivation
Compare Taylor’s Expansion of LHS and RHS
So for:
–
First expand the LHS
© Imperial College London
Derivation
• Now expand individual terms on RHS of our
expression:
• This gives:
© Imperial College London
Derivation
Equate both sides and solve for constants:
Similarly we can find C, D, E, F and C, D, E, F
© Imperial College London
Applying 6th Order Method
How to obtain results using derived method.
Produce a set of simultaneous equations and solve.
Find y’n+1 from old values and those just obtained.
Set x, yn , y’n etc. to new values.
Repeat above procedure with updated variables.
© Imperial College London
Circle Problem Example
Circle is produced for phase plane plot
© Imperial College London
Kepler’s Equations
• Solution obtained from
z-y plot reinforces fourth
order method results.
© Imperial College London
The End
But perhaps you want more?
Read our report
Visit our website:
http://www.doc.ic.ac.uk/~pb401/DE