:

## Stability Regions Of Runge-Kutta Solvers

Maple 13

This post gives equation of stability region of a Runge-Kutta ODE solver (Dormand-Prince45) as implicit domain equation and ODE-IVP.

Contents are:

• Runge-Kutta solver is implemented
• Solver is applied to test equation and growth factor is obtained
• Stability regions are defined by separating hyperplane (GF=growth factor=1)
• Stability regions are plotted by bruteforce approach (implicitplot)
• Boundary of stability regions are described as ODE and solved numerically
• ODE solutions are written to text files
• Output text files are rendered by TikZ

Stability Regions of Runge-Kutta Methods

 > restart;

Define Runge-Kutta Integrator

Dormand-Prince Butcher Tableau

 > DPA := [1/5, 3/10, 4/5, 8/9, 1, 1]: DPB := [[1/5, 0, 0, 0, 0, 0], [3/40, 9/40, 0, 0, 0, 0], [44/45, -56/15, 32/9, 0, 0, 0], [19372/6561, -25360/2187, 64448/6561, -212/729, 0, 0], [9017/3168, -355/33, 46732/5247, 49/176, -5103/18656, 0], [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]]:

Dormand-Prince Method

 > ODE45DP := proc(f,t,x,h)         description "integrates state derivative function 'f(t,x)' by time step 'h' using Dormand-Prince (RK) method";         local f0,f1,f2,f3,f4,f5,tnew,xnew;         f0 := f(t,x);         f1 := f(t+h*DPA[1], x+h*(DPB[1,1]*f0) );         f2 := f(t+h*DPA[2], x+h*(DPB[2,1]*f0+DPB[2,2]*f1) );         f3 := f(t+h*DPA[3], x+h*(DPB[3,1]*f0+DPB[3,2]*f1+DPB[3,3]*f2) );         f4 := f(t+h*DPA[4], x+h*(DPB[4,1]*f0+DPB[4,2]*f1+DPB[4,3]*f2+DPB[4,4]*f3) );         f5 := f(t+h*DPA[5], x+h*(DPB[5,1]*f0+DPB[5,2]*f1+DPB[5,3]*f2+DPB[5,4]*f3+DPB[5,5]*f4) );         xnew := x + h * ( DPB[6,1]*f0+DPB[6,2]*f1+DPB[6,3]*f2+DPB[6,4]*f3+DPB[6,5]*f4+DPB[6,6]*f5 );         tnew := t + h;                  (tnew,xnew); end proc:

Test Function and Growth Factor

Above method is applied to test problem  to determine growth factor for successive state updates

 > # define test function testfun := (t,x) -> lambda*x;
 (3.1)
 > # apply Dormand-Prince method to test function tnew,xnew := ODE45DP(testfun, t, x, h): xnew := simplify(xnew);
 (3.2)
 > # determine growth ratio growthratio := simplify(xnew/x); growthratio := algsubs(h*lambda=z,growthratio);
 (3.3)

check the order upto which growth ratio aggrees to actual growth ratio

 > taylor(exp(z),z,7);
 (3.4)
 > # express growth ratio in x,y coordinates growthratioxy := evalc(abs(subs(z=x+I*y,growthratio)));
 (3.5)

critical value of growth ratio is 1 below and above iteration converges or diverges, respectively

 > F := growthratioxy-1; F2 := growthratioxy^2-1;  # if square root is not wanted
 (3.6)

Plot growth ratio in as a graph (3D) and density+contour (2D) plot

 > plots[display]( {         plot3d( log(growthratioxy), x=-4..2, y=-4..4,style=patchnogrid),         plots[contourplot3d]( log(growthratioxy), x=-4..2, y=-4..4, contours=[0],grid=[64,64],color=black,thickness=4)},axes=boxed); plots[display]( {         plots[densityplot]( log(growthratioxy), x=-4..2, y=-4..4, grid=[64,64],axes=boxed, style=patchnogrid),              plots[implicitplot]( log(growthratioxy), x=-4..2, y=-4..4, grid=[64,64],axes=boxed, style=patchnogrid,scaling=unconstrained, thickness=3)} );

direction normal to  contour is obtained by gradient  so direction along contour is  .

 > G := VectorCalculus[Gradient](F,[x,y]): T := <-G[2],G[1]>/sqrt(G[1]^2+G[2]^2):

 > plots[display]({         plots[fieldplot]([G[1],G[2]], x=-4..4,y=-4..4,fieldstrength=log,axes=boxed,color=blue,grid=[20,20]),         plots[fieldplot]([T[1],T[2]], x=-4..4,y=-4..4,fieldstrength=log,axes=boxed,color=green,grid=[20,20]),         plots[implicitplot]( log(growthratioxy), x=-4..2, y=-4..4, grid=[64,64],axes=boxed, style=patchnogrid,scaling=constrained, thickness=3)});

Examining implicitplot figure by "point probe" tool, a convenient point on big loop is  with and inside point of  The small loops are located approximately at Loop starting positions are found via numerical solution, viz.

 > # determine x values at y=3.6 solve( subs(y=3.6,F2), x): smallBottomLoopX := select(type,map(evalf,[%]),realcons);
 (3.7)
 > # determine x values at y=-3.6 solve( subs(y=-3.6,F2), x): smallTopLoopX := select(type,map(evalf,[%]),realcons);
 (3.8)
 > min(smallTopLoopX); Statistics[Mean](smallTopLoopX);
 (3.9)

as expected rom figure, x values for  found similar

Define differential equation so as to move point  in direction  i.e.  and  also to detect closure of loop/trajectory augment states by

 > DEQx := diff(x(t),t)=subs(x=x(t),y=y(t),T[1]): DEQy := diff(y(t),t)=subs(x=x(t),y=y(t),T[2]): DEQa := diff(theta(t),t)=normal(diff( arctan(y(t)-cy,x(t)-cx), t));
 (3.10)
 >
 > StabilityRegion1 := dsolve({DEQx,DEQy,subs(cx=-2,cy=0,DEQa),x(0)=0,y(0)=0,theta(0)=0},'numeric',{x(t),y(t),theta(t)},events=[[theta(t)-2*Pi,halt]]); # StabilityRegion2 := dsolve({DEQx,DEQy,subs(cx=0.94,cy=-3.6,DEQa),x(0)=max(smallTopLoopX),y(0)=-3.6,theta(0)=0},'numeric',{x(t),y(t),theta(t)},events=[[theta(t)-2*Pi,halt]]); # StabilityRegion3 := dsolve({DEQx,DEQy,subs(cx=0.94,cy=+3.6,DEQa),x(0)=max(smallTopLoopX),y(0)=+3.6,theta(0)=0},'numeric',{x(t),y(t),theta(t)},events=[[theta(t)-2*Pi,halt]]); #
 (3.11)
 > StabilityRegion1(0.1); eval( [x(t),y(t),theta(t)], StabilityRegion1(0.1));
 (3.12)
 > region1plot := plots[odeplot](StabilityRegion1, [x(t),y(t)],t=0..17,numpoints=100,axes=boxed,scaling=constrained); region2plot := plots[odeplot](StabilityRegion2, [x(t),y(t)],t=0..2,numpoints=100,axes=boxed,scaling=constrained); region3plot := plots[odeplot](StabilityRegion3, [x(t),y(t)],t=0..2,numpoints=100,axes=boxed,scaling=constrained);
 (3.13)
 >
 > plots[display]({ region1plot, region2plot, region3plot });
 > writeLoopToFile := proc( fname, dsol, dt, tmax)         local T, fid, i, xy, val;         fid := fopen(fname,WRITE,TEXT);         T := 0;         dsol := eval(dsol); # do full evaluation         for i from 1 to round(tmax/dt) do                 val := dsol(T);                 xy := eval( [x(t),y(t)], eval(val));                 fprintf(fid, "%12.5f %12.5f\n", xy[1], xy[2]);                 T := T+dt;         end do;                 fclose(fid); end proc;
 (3.14)
 > writeLoopToFile( "Region1.dat", StabilityRegion1, 0.1, 16); writeLoopToFile( "Region2.dat", StabilityRegion2, 0.025, 1.4); writeLoopToFile( "Region3.dat", StabilityRegion3, 0.025, 1.4);

created loop coordinate data files are can be rendered in TikZ. An example rendering document is given below

\documentclass{article}

\usepackage{tikz}

\begin{document}

\begin{tikzpicture}

\draw[very thin,dashed,color=gray] (-3.9,-3.9) grid (3.9,3.9);

\draw[thick, red, ->] (-0.2,0) -- (4.2,0) node[right] {$\lambda_x$};

\draw[thick, red, ->] (0,-1.2) -- (0,4.2) node[above] {$\lambda_y$};

\draw[green] node at (0., 5.) {\textsc{Stability regions for Dormand-Prince 45 method}};

\draw[thick, blue, fill=blue,fill opacity=0.5] plot[smooth] file{Region1.dat};

\draw[thick, blue, fill=blue,fill opacity=0.5] plot[smooth] file{Region2.dat};

\draw[thick, blue, fill=blue,fill opacity=0.5] plot[smooth] file{Region3.dat};

\end{tikzpicture}

\end{document}

Result of rendering (using MikTex) is given below

 >
 >
 >