Overview to adaptive shooting methods for boundary value problem

Implementation of numerical methods related to nonlinear boundary value problems of ordinary differential equations

Josef Otta

University of West Bohemia, Czech Republic

Overview

Shooting method for ODE

1D method

Introduction

Algorithm

Numerical examples and applications

2D method

Principle of method

Algorithm and principle of adaptivity

Numerical examples

Conclusion

Shooting Method for boundary value problem of ODE

Principle of the method

One of the classic numerical method for solving non-linear boundary value problem for ordinary differential equation is Shooting method -  a method transforming problem of solution boundary value problem to solving sequence of solutions of ODE with initial contitions. Let's describe it on simple example of Dirichlet problem for ODE of 2nd order.

u''=f(t, u, u'),t ∈ ⟨a, 
          b⟩
u(a) =(Subscript[η, a])
u(b) =(Subscript[η, b])

where f is continuous function.

If the function f is Lipschitz continuous than Cauchy problem has unique solution with any initial value at t=a.

u''=f(t, u, u'),t ∈ ⟨a, 
          b⟩
\(u \((a)\) = η\)
u' (a) = ξ

Our goal is to find such values of  η and ξ  that solution of initial value problem is also solution of boundary value problem. We can transform this problem to problem of finding roots of non-linear functional F(u)=0. This functional  represents behaviour of computed solution with respect to boundary conditions.

1D shooting method

Principle of algorithm

Usage of adaptivity

To find the roots of functional F we have to know something about its behaviour. Since the functional F is studied by its numerical values we can use built-in function Plot which is primarily used for drawing graphs of functions with 1 independent variable. Graph of the function is approximated by graph of piecewise-linear function where Mathematica uses following adaptive procedure.
	First of all, equidistant points are computed. Number of these points is determined by parameter PlotPoints. Then sample points are filled in places where succesive points contain higher angle than  parameter MaxBend. Recursion limit of this procedure is determined by parameter PlotDivision.
	In sysem Mathematica is everything represented by abstract data-structure expression. Thanks to this we can extract coordinates of sampling points from plot function and used them for future computation. For illustration how subdivision method works is  picture below presented.

graf=Plot[Sin[x]/x,{x, -10, 10},PlotPoints -> 5,PlotDivision -> 10,DisplayFunction -> Identity];
body = 
      Cases[graf, _Line, ∞]; 
pts = 
      ListPlot[Partition[Cases[body, _?NumberQ, ∞], 2], 
        PlotJoined -> False, 
        PlotStyle -> {PointSize[0.020], Hue[0.4,  0.4,  0.6]}, 
        DisplayFunction -> Identity]; 
Show[graf, 
      pts, DisplayFunction -> $DisplayFunction];

1D shooting method

Principle of algorithm

Steps of algorithm

1. Definition of the equation and parameters

2. Define the functional f(u)

3. Extract the coordinates of points from plotted functional

4. Find the intervals where the sign of values is changing

5. Compute roots with Brendt's algorithm

6. Solve equation with initial values computed initial conditions

7. Show solution and compute error norm

1D shooting method

Numerical example

This method was developed for solving piecewise linear differential equations of second order with Dirichlet, Neumann and Robin conditions.

\(u'' \((t)\) - \(λ\_+\)\ \(Max \((u \((t)\), 
                  0)\)\) + \(λ\_-\)\ \(Max \((\(-\(u \((t)\)\)\), 
                  0)\)\) + Sin \((u \((t)\))\)\  = \ 0, \ 
        t ∈ ⟨0, π⟩\)
(Subscript[u, a]) =α u (a) + β u' (a);    |α|\(+\ \(\(|\)\(β\)\(|\)\(\ \)\(\(≠\)\(0\)\)\)\)
(Subscript[u, b]) =γ u (b) + δ u' (b);    |γ|\(+\ \(\(|\)\(δ\)\(|\)\(\ \)\(\(≠\)\(0\)\)\)\)

Since the equation is nonlinear we can expect more solutions. For that reason we are using our adaptive shooting method.

We define functional F(u) as

F (u (η)):=γ(Subscript[u, a])(b, η)+δ u' (t, η)-(Subscript[u, b])

1D shooting method

Numerical example

Dirichlet problem of nonlinear equation with jumping nonlinearities

\(u'' \((t)\) - \(λ\_+\)\ \(Max(u(t), 
              0)\) + \(λ\_-\)\ \(Max(\(-\(u(t)\)\), 0)\) + 
          Sin(u(t))\  = \ 0, \ 
      t ∈ ⟨0, π⟩\)
u(0) = 0
u(π)=0

\(where\ \(λ\_+\) = \(1\ and\ \(λ\_-\) = 16. \)\)

({{PlotPoints, 150}, {"PlotDivision", 10}, {Sampling points, 261}, {"Number of roots:", \(\*"""" 5\)}, {"Timing", ""}, {"Plot", 1.297` Second}, {"FindRoot", 0.18800000000000017` Second}, {"NDSolve", 0.03100000000000014` Second}, {"Total", 1.5160000000000002` Second}, {"Error", ""}, {"Maximal Error:", 3.2180704993199256`*^-11}})

1D shooting method

Application of this method

Approximation of Fučík spectrum

\(u'' \((t)\) - \(λ\_+\)\ \(Max(u(t), 
              0)\) + \(λ\_-\)\ \(Max(\(-\(u(t)\)\), 0)\)\  = \ 0, \ 
      t ∈ ⟨0, π⟩\)
u(0) = 0
u(π)=0

\(where\ [\(λ\_+\), \(λ\_-\)] ∈ \ ℝ\^2\)

Now, we use ContourPlot which shows the contours of the  function f : (R^2)⟶(R^1)where independent variables are \(λ\_+\)and \(λ\_-\), \(dependent\ variable\ is\ \ value\ of\ the\ \(\(functional\)\(.\)\)\)

We can obtain Fučík's spectrum by drawing zero contours of our function.

2D shooting method

Introduction

This method was developed for solving piecewise-linear differential equations of second order with periodic boundary conditions. Our goal was to find all existing solutions of following equation

u''(t)+β u'(t)+\(λ\ \(\(u\^+\) \((t)\)\)\)+u (t)=1 + ϵ cos (2 t)
u (t) = u (t + T),   T > 
          0.

For solution of this problem we used also shooting method but with 2 parameters.

u''(t)+β u'(t)+\(λ\ \(\(u\^+\) \((t)\)\)\)+u (t)=1 + ϵ cos (2 t)
u (0) = ξ, u' (0) = η,

We define functional F(u) as

F (u (ξ, η)):=\(\@\(\((u' \((T)\) - u' \((0)\))\)\^2 + \((u(T) - u(0))\)\^2\)\)

2D shooting method

Principle of algorithm (1/2)

Usage of adaptivity

Purpose of this method is computing all existing global minimas of this nonnegative function. If we are more exact we are trying to find points where functional has zero values. We are making adaptive sampling in the areas where we expect minimum of the functional. So we do not make fine sampling in areas where the information would be redundant. Illustrative scheme is shown at the following picture. The principle of this method is to make a system of subareas where a rule of inclusion child ⊂ parent is valid.

2D shooting method

Principle of algorithm

Basic idea of finding all of the minimas (1/2)

Let's imagine a simple example of a function of 2 variables.

f(x, y)=sin(2 (Sqrt[(x^2) + (y^2)])) (x^2) + (y^2) + 4

f[x_, y_] := 
      x^2*Sin[2 √((x)^2 + ( y)^2)] + (y^2) + 
        4
Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, PlotPoints -> 30, 
      Mesh -> False, 
      ColorFunction -> (Hue[(#^1.4), 0.6, 1] &)];

2D shooting method

Principle of algorithm

Basic idea of finding all the minimas (2/2)

cont = 
      ContourPlot[f[x, y], {x, -2, 2}, {y, -2, 2}, 
        ColorFunction -> (Hue[1, 1 - #, 1 - 0.7  #^4] &), 
        PlotPoints -> 50, Contours -> 10, 
        ContourLines -> True];

Contour  Plot

Choosing the areas of interests

Enclosing areas to rectangles

2D shooting method

Convergence criterion

Lemma:  Sufficient condition of convergence of quntile sequence
	Sequnce of quantiles  {q(n)}_(n = 1)^∞ is convergent if  condition

q(n) ≤ ((q(n - 1))/M),    M > 
      1

is satisfied.

High estimate of subareas size

▫	computed areas have diagonal character
▫	subareas are enclosed to rectangles -- the speed of reducing size of subareas is slower.
▫	high estimate of tolerance value after n steps is (((((q(1))/M)))^n)
▫	size of reduced areas  -- size of original twodimensional interval

High estimate of the speed  of contraction  of subareas  size

▫	subareas are approaching to rectangles.
▫	size of subareas after n steps |int| ((((k/100)))^n)
▫	high estimate of tolerance value after n steps (((((q(1))/M)))^n)

2D shooting method

Principle of algorithm

Steps of algorithm

1. Definition  of q(0), M, k,n,l, tol
2. i:=1
3. Equidistant sampling of the points of functional F
4. q(i) := k-percent quantile from sampled points
5. If q(i)>=((q (i - 1))/M)  then  q(i) := q (i - 1)/M
6. Clustering the points into subareas Pij ; m = number of subareas;
7. For j = 1..m
	(a) If Pij  is not a point, 
		then
			i. Enclosing areas to rectangles
			ii. saving quantile q(i)
		 else save as initial point for function FindMinmum;
	(b) number of all the subareas := pij ;
	(c) If pij > 0, then  i:=i+1; for all areas   goto 2,
      else goto 3;
8. OBL := {Pn,i | subarea sampled in n-th recursive step}      
9. find minimas at set of areas OBL and set {Pbij} in high precision
10. If computed minimas (Subscript[m, i]) < tol, then (Subscript[m, i]) is solution

2D shooting method

Finding the minimas

After adaptive sampling method of the area  int we get the subareas. Now we suppose that they are small enough and that the solutions differ from themselves more than the diameter of the computed subareas.

If we suppose this fact, we can use the methods for finding local minimum on these subareas. We are computing with FindMinimum function in high precision. Then all of the computed local minimas are considered to be global ones.

2D shooting method

Samle output

cas = 
      Timing[Shootingmethod2D[{64, 64, 1}, 2, 50, 
          0.15, {{-2, 2}, {-2, 2}}, 1, 10]]

2D shooting method

Numerical approximation of solution

We can apply this method to solving  boundary value problem for ordinary differential equation of the second order with jumping non-linearities describing onedimensional model of suspended bridge in form

u''(t)+β u'(t)+\(50 \(\( u\^+\)(t)\)\)+u(t)=2 cos(2 t) + 1
u(t) = u(t + 4  π)

Solution with initial condition u(0)=0.1698 a u'(0)=0.4556.

Solution with initial condition u(0) = 0.1681  and u'(0) =0.5035

Solution with initial conditions u(0) = 0.1534 and  u'(0) =0.8242

Conclusion

Proposal of adaptive 1D method

Adaptive 2D method

could help to answer some open quesions concerning bifurcation problems

Universal  sampling

\(which\ could\ guarantee\ that\ all\ of\ existing\ minimas\ can\ be\ found\ \((for\ function\ from\ C\^2\ \(\(space\)\(.\)\))\)\)

This work was made within project Realization of interactive-informational portal for scientific-technical apllications, number 1N04078, financed by Ministry of Education of The Czech Republic and solved in University of West Bohemia in Plzeň, Czech Republic.

References

[1] Coddington,E.A.,Levinson,N.,Theory of Ordinary Diferential Equations.McGraw–Hi ll,New York 1955.

[2] Čepička,J.,Numerické experimenty v nelineárních úlohách.Disertační práce,Fakulta aplikovaných věd,Západočeská univerzita v Plzni,Plzeň 2001.

[3] Doole,S.H.,Hogan,S.J.,A Piecewise Linear Suspension Bridge Model:Nonlinear Dynamics and Orbit Continuation.In proceedings:Dynamics and Stablity of Systems 11(1),1996.

[4] Míka,S.,Přikryl,P.Numerické metody řešení obyčejných diferenciálních rovnic,Vydavatelství ZČU,Plzeň 1994.

[5] Otta,J.,Numerické řešení nelineárních okrajových úloh metodou střelby v systému Mathematica. In proceedings: Proceedings od the 4th International Conference APLIMAT 2005,Bratislava,Slovenská republika,2005,s.579-584. ISBN 80-969264-3-8.

[6] Přikryl,P.,Brandner,M.,Numerické metody II,Vydavatelství ZČU,Plzeň 2000.

[7] Rychlík,J.,Programovací techniky,Nakladatelství Kopp,České Budějovice 1994.

[8] Wolfram,S.,Mathematica 5.0 – Advanced Documentation,Wolfram Resarch,Inc.,2003.

[9] Wolfram,S.,Mathematica Book 4th ed.,Wolfram Media/Cambridge University Press,1999.

[10] Žára,J.,Beneš,B.,Felkel,P.,Moderní počítačová grafika,Computer Press,Praha 1998.