User:Ramiamro

=Project Report for User:Ramiamro=

For Introduction to Numerical Analysis, Fall 2010.

Introduction
the project intended to clarify the concept of the stability of ordinary differential equations' Methods, and especially Runge-Kutta methods of order four RK4,it is also discussing the shape of the absulute stability regions for each method, many examples, quizes, and exercises have been established on wikiversity and wikipediastiff equation, this establishment focused on users who have little experience with stability, such that thay can grasb the concept of absolute stability as well as A/L- stablility, more over the project studies particular examples of finding the characteristic polynomial for the multistep methods.

To operate numerical methods on stiff equations and get good approximation for the exact solution; methods with special stability are required, like A-stable or L-stable methods, if the method is applied with h>0 to the test equation then for A-stable method, it will produce approximation function $$R(\lambda h)$$ to the exact solution which will be proportional to $$ exp(\lambda h)$$, the approximation function will be bounded by 1 in the left half of the complex plane, which gave us the stability region, and so we can chose the step size which suitable for a certain kind of equation using that method.

Initial Experience

 * During the course and in courses before I introduced to many types of stabilities, but each time I have those we talked about similar diffenetions of stability in different topics, like stability for ODE's and stability of nonlinear systems, as well the stability for physical systems.
 * To study stability it is important to define the area you study in, the stability diffinitions are similar but at the same time time they are significantly different in details; in numerical analysis we talk about the stability of the method when we apply it to a certain type of equations. i.e. the test equation, our cocentration during this course was on the stiff equation, which is a differential equation that shows wide variance in the results as a result of small changes in the argument (the independent variable), some of the numerical methods shows very good stability when treating such kind of equations, specially the implicit methods, we will see during the project how the explicit methods Rk4, and others how they behave when they employed to simulate the exact solutions of those equations.these modifications manifisted in posting some diffinitions and an example on the wikipedia page stiff equation, and some in quizes and excercices on my wikiversity user  page, these changes took place because I found they are not clearly discussed on the wikipedia page, although my changes my be in the hand or the minds of others, but I beleive I made a better way of  understandig the concept of stability for the Numerical methods more clearer, and easily comprehended especially by using the conceptual quiz I posted on wikiversity down.

Motivation
During the couse, we covered so many topics, and I barely remeber that any class passed without asking or discussing the error propagations in the numerical method which is implicitly implies whether the numerical method is stable or not, on the other hand during the second half of the course these cries increased and I felt that the stability topic is not clear for some of the students, due to the time limitations of the class and the poorley expressed ideas on wikipedia I made my point to clarify this topic more and make it easier to comprehend, so I established my project with a hard base to study the stability of the numerical method such that I introduced the reader to basic diffinitions in this topic, so he can cope up with the flow of the material later on, and made sure the reader will not losely use the concept later on, then moved to show some examples "mine and others on wikipedia" so that he/she get experience with the difinitions

Proposed Changes
my proposed changes were considering,changes to the following subtopics:
 * Stability regions for RK methods.
 * Stability regions for multistep method.
 * The characteristic functions and how to find it.
 * Solving a system of differential equations using Rk4 method.
 * plotting the stability regions.
 * Quizes related to the topic: stability of ODEs' methods.
 * Excercises related to the same topic.

Actual Changes
I almost covered my plan, but I did too much in RK methods than the implicit multistep analysis, one thing that I didn't cover is the application of RK4 to a system of differential equations, the idea is fancy to be discussed on wikipedia, since very little details written about it, and it is worthy to show some hidden details, most books they just give the general formula to apply that method to a system of differential equations.

I starrted my project with a number of definitions such as consistency, convergense, and absolute stability,then followed by a theorem which uses those concepts, the theorem shows some results concerning stability and convergence of the numerical methods both single and multistep methods, this is an excution for what I proposed about the multistep method, in addition to excercises that covers finding the characteristic polynomial for RK2 method, and a nother method.

The rest of my proposed changes has been covered and an example to show how to find the stability polynomial for five types of RK4, including the classical case, A matlab code was shown and tested to generate the stability region for the five cases,In fact I generated the stability regions for the five cases and compared the with other resources, An image of the generated stability region was published on stiff equation, all the changes have been published to the wikipedia. concerning the quizes part I managed to write 13 conceptual quizes, since finding numerical quiz in this topic is really long way, so I just posted the 13 questions on wikiversity below, any user can also access them from user:Mjmohio, under Numerical ordinary differential equations

Excercises have been developed, I published basically tow questions, divided them to five guiding excercises, and arranged them in chronological order in terms of solving for the stability region, and testing wether the the numerical method is stable or not.During the exercises I clarified the L-stability concept by solving and discussing an L-stable method,

Definitions
let $$\tau_i$$ is the local truncation error, and k is the number of time steps, then:   The numerical method is consistent with a differential equation if


 * $$\lim_{h \to 0} \operatorname{max}|\tau_i|=0$$ over $$1\leqslant i \leqslant k$$.

according to this definition Euler's method is consistent.

 A numerical method is said be convergent with respect to differential equation if


 * $$\lim_{h \to 0}|x(t_i)-y_i|=0$$ over $$1\leqslant i \leqslant k$$;

where $$y_i$$ is the approximation for $$x(t_i)$$.

 A numerical method is stable if small change in the initial conditions or data, produce a correspondingly small change in the subsequent approximations. 

Theorem: For an initial value problem
 * $$x'=f(t,x)$$ with $$t\in[t_0,t_0+\alpha]$$

and certain initial conditions, $$(t_0,x_0)$$, let us consider a numerical method of the form
 * $$(y_0=x_0)$$ and $$y_{i+1}=y_i+h\phi(t_i,y_i,h)$$.

If there exists a value $$h>0$$ such that it is continuous on the iterative domain, $$\Omega$$ and if there exists an $$L>0$$ such that
 * $$|\phi(t,y,h)-\phi(t,y^*,h)|\leqslant L |y-y^*|$$ for all $$(t,y,h),\,(t,y^*,h)\in \Omega$$,

the method fulfills the Lipschitz condition, and it is stable and convergent if and only if it is consistent. That is,
 * $$\phi(t,x,0)=f(t,x)$$ for all $$t \in \Omega$$.

For a similar argument, one can deduce the following for multi- step methods:   The method is stable if and only if all roots, $$\lambda$$, of the characteristic polynomial satisfy
 * $$|\lambda|\leqslant 1$$,

and any root of
 * $$|\lambda|=1$$

is simple root.  one more result is that if the method is consistent with the differential equation, the method is stable if and only if it it is convergent.  see

stability polynomials of Runge-Kutta methods
The Runge–Kutta methods are very usefull in solving systems of differential equations, it has wide applications for the scientists and the engineers, as well as for the economical models, the recognized with their practical accuracy where we can use and get very good results and approximations when solving an ODE problem, RK has the general form :$$ y_{n+1} = y_n + h\sum_{i=1}^s b_i k_i$$,

Example:finding the stability polynomial for RK4's methods
for RK4" case$$2_a$$", which characterizes by $$c_2= \frac {-1}{4}$$ which has the form: the stability region is found by applying the method to the linear test equation $$y'=\lambda y$$
 * $$ y_{n+1} = y_n + \frac{h}{6}k_1+0k_2+\frac {2k_3}{3}+\frac {k_4}{6} $$
 * the stability region is found by applying the method to the linear test equation $$y'=\lambda y$$


 * $$ \displaystyle \ k_1=hf(t_n,y_n)$$
 * $$k_2=hf(t_n+\frac{h}{2},y_n-\frac{k_1}{2})$$
 * $$k_3=hf(t_n+\frac{h}{2},y_n+\frac{3 k_1}{4}-\frac{ k_2}{4})$$
 * $$k_4=hf(t_n+h,y_n- 2k_1+k_2+2k_3)$$

using the linearized equation $$f(t,y)=\lambda y$$, we get
 * $$k_1=h\lambda y_n$$
 * $$k_2=h\lambda(1-\frac{h\lambda}{2})y_n$$
 * $$k_3=h\lambda(1+\frac{3 h\lambda}{4}-\frac { h \lambda}{4}(1-\frac{h\lambda}{2}))y_n$$
 * $$k_4=h\lambda(1-2h\lambda+ h\lambda(1-\frac{h\lambda}{2})+2h\lambda(1+\frac{3 h\lambda}{4}-\frac { h \lambda}{4}(1-\frac{h\lambda}{2})))))y_n$$

substitue these back in $$y_{n+1}$$, yields $$y_{n+1} = [1+(h\lambda)+\frac{(h\lambda)^2}{2}+\frac{(h\lambda)^3}{6}+\frac{(h\lambda)^4}{24}]y_n=R(h\lambda)y_n $$
 * and so the characteristic polynomial

$$ P(z)=z-R(z); z=h\lambda $$ for the absolute stability region for this method, set |R(z)|<1, and so we get the region in figue1. case$$2_a$$

the table below shows the final forms for the stability function for different forms of RK4, these RK4's are different in the values of $$ b_j$$, and they are fullfilling the consistencty requirement for the method i.e :$$\sum_{j=1}^{i-1} a_{ij} = c_i\ \mathrm{for}\ i=2, \ldots, s.$$

Plotting the stability region
In order to plot the stability region, we can set the stability function to be bounded by 1 and solve for the values of z, then draw z in the complex plane. Since R(z) is the unit circle in the complex plane, each point on the boundary can be represented as $$e^{i \theta}$$ and so by changing $$\theta$$ over the interval$$[0,2\pi]$$, we can draw the boundaries of that region. The following OCTAVE/Matlab code does this by plotting contour curves until reaching the boudaries: The figure at right shows the absolute stability regions for RK4 cases which is tabulated above

Quizes
{ the A-stability is characterized by: - all the points in the left half of the complex plane. - it a special case of the absolute stability. - inludes all $$\{ z \in \mathbb{C} | \mathrm{Re}(z) < 0 \}$$, + all of the above
 * type=""}

{Euler method is: - A-stable + not A-stable
 * type=""}

{the absolute stability region for RK4 is greater than the A-stability region for the same method:
 * type=""}

+ false - True - True for certain values of $$ b_j$$.

{A numerical method is stable if : + small change in the initial condition will produce small change in subsequent steps. - small change in the initial condition will produce huge change in subsequent steps. - big change in the initial conditions produce oscillatory solution at the end.
 * type=""}

{$$ \lambda =1$$ is always a root of the characteristic polynomial of the multistep method:
 * type=""}

+ True - False - usually not the case

{the root of the C.P. can be real or complex,and the method still be stable: + True - False
 * type=""}

{A multi- step method is strongly stable if : - $$ \lambda =1$$ is the only root of magnitude 1. - all other roots has magnitude <1 + the first and the second sentence - just one root inside the unit circle if it has more than one.
 * type=""}

{if more than one root has magnitude equal to 1, and the others are less than one, the method is - strongly stable. - A-stable + weakly stable. - Unstable.
 * type=""}

{The numerical method is unstable if - $$ \lambda > 1$$ for at least one root - $$ Re(\lambda) > 1$$ for at least one root - $$ Re(\lambda) \leqslant 1$$ for at least one root + the firs and the second sentences.
 * type=""}

{the absolute stabity region for the explicit Euler's method is + unit cicle in the complex plane, its center shifted to the left, by one unit. - unit cicle in the complex plane, its center shifted to the rigt, by one unit. - the whole left side of the complex plane. - the method is unstable.
 * type=""}

{ None of the RK methods is A-Stable: + True - False
 * type=""}

{All explicit methods are : - A-Stable + Not A-Stable. - Unstable. - None of the above
 * type=""}

{implicit multistep methods are A-stable if the have order at most -1 +2 -3 -they are always A-stable.
 * type=""}

Ex:1
find the stability function for RK2 which is given by: $$ y_{n+1}=y_{n}+hf(t+\frac {h}{2},y_{n}+\frac {hf(t_n,y_n)}{2})$$

Solution: applying this method to the test equation $$ y'=\lambda y$$
 * we get $$ \begin{align}

y_{n+1} &=y_n+h \lambda y_n+\frac {(h \lambda)^2}{2}y_n \\ \Rightarrow \end{align}$$ the stability polynomial $$ \phi(z)=1+h \lambda +\frac {(h \lambda)^2}{2}$$

Ex:2
find the absolute stability region for RK2. Solution: by setting $$ z = h \lambda $$
 * $$ \Rightarrow $$ the abs.stability region is given by $$\{z\in \mathbb{C}||1+ h \lambda + \frac {(h \lambda)^2}{2}| < 1 \}$$

Ex:3
find the characteristic polynomial for RK2. Solution: it is $$ r^{i+1}=(1+z+\frac {z^2}{2})r^i$$ divide both sides of the equation by $$ r^i;r \ne 0$$ you get $$ r^{i+1}=(1+z+\frac {z^2}{2})r^i$$
 * $$ \Rightarrow \phi(r,z)=r-1+z+\frac {z^2}{2}$$

Ex:4
is RK2 stable, if it is what type of stability. Solution: you get $$ r^{i+1}=(1+z+\frac {z^2}{2})r^i$$ by setting z=0,
 * $$ \Rightarrow r=1$$
 * so the method is strongly stable since r=1, is the only root, and has a value of 1.

Ex:5
Determine the stability of Back ward Euler method. Solution: $$ y_{n+1}=y_n+hf(t_{n+1},y_{n+1})$$
 * by applying this method to the test equation
 * then $$ y_{n+1}=y_n+h\lambda y_{n+1}$$
 * and so $$ y_{n+1}=\frac{1}{1+h\lambda }y_{n}$$
 * call $$ z=h\lambda$$ and so:
 * $$ G(z)\longrightarrow 0$$ as $$ Re(z)\longrightarrow \infty$$
 * and so this is L-Stable method when applied to stiff equation.

Conclusions
During this project, I searched about five resources and other internet websites to weave my ideas and make them as clear as possible, I used examples to clarify the definitions and the theory, I found that Rk methods of order one and the higher orders are cositent with stiff problems, and so they are convergent, see the proof on , and so when we use RK methods we will not be worried if the method is consistent or not. Using any method to simulate the exact solution for any kind of ODE's is dependent on the absolute stability region, and so the choice of the time step will depend on the roots of the characteristic polynomial, if the complex number $$ z=h\lambda$$ should satisfy the condition that the absolut value of the stability polynomial should be strictly less than one, otherwise the method will not be stable, and the results will be far from the accuracy.

On the other hand, theory should be clarified by examples, it wasn't easy to find the stability polynomial for Rk methods, one should be very carefull when use the parantheses, see the example above, missing of one parantheses will leed to a different stability polynomial and so leads to a different stability region, which in role means different method.

Conceptual quiz was very neat to clarify the concepts, and so I included so many of them, they expressing the definition of the absolute stability, A-stability, strong and weak stability, and L-stability; see exercises above, they inclded more about stability properties of RK methods, and recruiting the attention of the reader to some general rules and corollaries that is not of first look to catch and comprehend. I finished my main modifications on wikipedia by showing a sample matlab code how to plot the stability regions,the code is working efficiently, and gives accurate results, see.

Note: it is worthy to mention, that during my writings I noticed that some resources define the stability region to be $$\{z \in \mathbb{C} | |\phi(z)|\leqslant 1 \}$$. But on [w:stiff equation] the definition doesn't include 1