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):

plot gradient (blue) and gradient normal (green) vector fields

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);

Warning, cannot evaluate the solution further right of 15.587999, event #1 triggered a halt

 

 

Warning, cannot evaluate the solution further right of 1.3594372, event #1 triggered a halt

 

 

Warning, cannot evaluate the solution further right of 1.3594361, event #1 triggered a halt

 

(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);

Warning, cannot evaluate the solution further right of 15.587999, event #1 triggered a halt

 

Warning, cannot evaluate the solution further right of 15.587999, event #1 triggered a halt

 

Warning, cannot evaluate the solution further right of 15.587999, event #1 triggered a halt

 

Warning, cannot evaluate the solution further right of 15.587999, event #1 triggered a halt

 

Warning, cannot evaluate the solution further right of 1.3594372, event #1 triggered a halt

 

Warning, cannot evaluate the solution further right of 1.3594361, event #1 triggered a halt

 

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

 



Download ODE45_ev.mw

plotStabRegionInTikZ.pdf

 


Please Wait...