Parallel Computing Toolkit and Numerical Experiments

Petr Girg
University of West Bohemia, Pilsen, Czech Republic

Special thanks

Jan Čepička, pioneer in Numerical Experiments ;
Jan Frisch,     system designer and administrator of  grid -   and
                                  distributed - 
          computing resources in our department.

Introduction

In our talk  we would like to briefly acquaint the audience with prototype problems studied in our workgroup, to show how to use Mathematica in our problems, and finally we would like to suggest improvements in Mathematica which would make our work more effective. 

Numerous phenomena in physics, chemistry, biology, and economy can be modelled by boundary-value problems for differential equations. These boundary-value problems (BVPs for short) can often be written in an operator form

L (u) + H (u,  p) = f ,

here

L (linear or nonlinear) mapping between appropriate function spaces X and Y ;

H: X ×  (ℝ^n)→ Y being an nonlinear operator ;

f ∈ Y being fixed ;

p ∈ (ℝ^n) being parameters . 

Solvability, multiplicity and bifurcation of solutions are of particular interest for people from praxis. Note that the BVPs depending on parameters describe such important phenomena as resonance, in which the real systems usually undergo drastic changes leading to their collapse (suitable small perturbations in external forcings effect in big changes in their corresponding responses).

Examples of problems:

Mathematical models of stationary processes from nonlinear heat trasfer,
nonlinear reaction-diffusion theory, magnetohydrodynamics and glaceology :

- div(\(TraditionalForm\`\(\(|\)\(grad\ u\)\( | \^\(p - 
              2\)\)\(\ \)\(grad\)\(\ \)\)\)u) - λ \(TraditionalForm\`\(\(|\)\(u\)\( | \^\(p - 2\)\)\)\)u + g(λ, x, u)= f ,          in Ω ;
                                                             u=0,            on ∂Ω 
or its one-dimensional variant

- (|u'((|)^(p- 2))u')' - λ \(TraditionalForm\`\(\(|\)\(u\)\( | \^\(p - 2\)\)\)\)u + g(λ, x, u)= f,               in  (0,1);
                     u(0)=0,  u(1)=0 . 

Mathematical models of one-degree-freedom oscillator:

u'' + λ u + g(λ, x, u, u' ) = f

subject to the preriodic b.c. :

u(0)=u(π),  u'(0)=u'(π)

or subject to the Dirichlet b.c. :

u(0)=0 , u(π)=0

or the Neuman b.c.:

u'(0)=0 ,  u'(π)=0 .

With more general nonlinearities involved in the equations, one looses idea what can be proved or disproved. For that reason, numerical experiments become an indispensable tool of nonlinear analyst nowadays in the following way. In order to get at least rough idea about the qualitative behavior, one performs a set of numerical experiments by discretizing over the domain of parameters. This approach leads to an enormous amount of BVPs usually. To solve such a number of BVPs would hardly be possible without parallelization. The parallelization with respect to the discretized set of parameters can be considered as a pure data-parallel approach (same code, different data). Parallel Computing Toolkit (PCT for short) brings data-parallel methods into Mathematica efficiently.

Example of an numerical experiment:

Given 1<p, let us investigate the set of all (λ, u) ∈ "\[DoubleStruckCapitalR]"×(W^(1,p))(0, 1) satisfying the following

- (|u'((|)^(p- 2))u')' - λ \(TraditionalForm\`\(\(|\)\(u\)\( | \^\(p - 2\)\)\)\)u = (f^T)+ a (Subscript[φ, 1]) ,               in  (0,1);
                     u(0)=0,  u(1)=0 . 
                     
boundary value problem with λ being close to (Subscript[λ, 1]).

Here (Subscript[λ, 1]) ∈ "\[DoubleStruckCapitalR]" and (Subscript[φ, 1]) ∈ (W^(1,p))(0, 1) satisfy

- (| (Subscript[φ, 1])'((|)^(p- 2)) (Subscript[φ, 1])')' - (Subscript[λ, 1]) | (Subscript[φ, 1])\( | \^\(p - 2\)\) (Subscript[φ, 1]) = 0 ,               in  (0,1);
                      (Subscript[φ, 1])(0)=0,   (Subscript[φ, 1])(1)=0 ;
             (∫_0^1)  | (Subscript[φ, 1])\( | \^p\)=1.
Note that this problem is related with the so-called nonlinear Fredholm alternative.
Fredholm alternative is the necessary and sufficient condition for the solvability of 
abstract equations with compact linear operators. In 1960-ties atempts to extend 
this theory to nonlinear but still homogeneous operators have been started. \(TraditionalForm\`\*"""\)It is easy 
to see that the differential operator u ↦ -(|u'((|)^(p- 2))u')' - λ \(TraditionalForm\`\(\(|\)\(u\)\( | \^\(p - 2\)\)\)\)u  \(TraditionalForm\`\*"""\)is (p-1)-homogeneous:              
     
-(|(a u)'((|)^(p- 2))(a u)')' - λ \(\(|\)\((a\ u)\)\( | \^\(p - 2\)\)\((a\ u)\)\) =-\(TraditionalForm\`\(\(|\)\(a\)\( | \^\(p - 2\)\)\(a\)\(\ \)\)\){(|u'((|)^(p- 2))u')' - λ \(TraditionalForm\`\(\(|\)\(u\)\( | \^\(p - 2\)\)\)\)u}




Parallelization with NDSolve

For mathematically correct setting of the problem and thorough discussion of numerical experiments, we kindly invite the reader to see [3], [4]. 

It is well known fact that the shooting method is very robust approach of solving BVPs for ODE's. Corresponding initial-value problem for the first-order system of ODE's:

u'  =  |v((|)^(p'-2))v ;
v'  = -λ |u\( | \^\(p - 2\)\)u  + (f^T)+ a (Subscript[φ, 1]) ;
(Subscript[φ, 1])'=  |w((|)^(p'-2))w ;
w' =  λ |(Subscript[φ, 1])((|)^(p - 2))(Subscript[φ, 1])  

subject to

u(0)  = 0
v(0)  =  |u(0)\( | \^\(p - 2\)\)u(0)  = |α\( | \^\(p - 2\)\)α   
(Subscript[φ, 1])(0) = 0
w(0)  = (((Subscript[π, p]) )^(p - 1))

Resulting system contains nonlinearities which are not all Lipschitz continuous. For that reason, uniqueness of solutions is not guarranteed. We speak of numerical experiment rather than of numerical computation of solutions.

We are ready to apply numerical methods now.

General philosophy of problem solving in Mathematica

Use build-in Mathematica functions as much as possible. 

General philosophy of data-parallel approach:

Divide the problem into independent parts where different data are  processed by the same code. 

Solution:

Let u(λ, α, a, 1) denote the value of the solution u to the initial boudary value problem.

It is not difficult to define corresponding numerical function in Mathematica using 
NDSolve . 


We parallelize with respect to the parameter a ∈ "\[DoubleStruckCapitalR]", i.e., we consider set of discrete values
 (Subscript[a, i]) ∈ "\[DoubleStruckCapitalR]" ,   i ∈ {1, 2, 3, ..., n}, where n ∈ ℕ is an multiple of the number of processors available, typically. 


Value of (Subscript[a, i]) given, we use ContourPlot to function  u(#1, #2, (Subscript[a, i]), 1)&  with an option Contours -> {0}.
Complexity of this step is O((PlotPoints^2)) passed to ContourPlot. 


Nodes of contour lines are then refined using FindRoot ;  complexity being O(PlotPoints).

See figures attached in postscript files Fredholm.ps, LandesmanLazer.ps, Oscilatory.ps.




Parallelization for symbolic code

Given 1<p, let us find (λ, u) ∈ "\[DoubleStruckCapitalR]"×(W^(1,p))(0, 1) satisfying the following

- (|u'((|)^(p- 2))u')' - λ \(TraditionalForm\`\(\(|\)\(u\)\( | \^\(p - 2\)\)\)\)u = (f^T)+ a (Subscript[φ, 1]) ,               in  (0,1);
                     u(0)=0,  u(1)=0 . 
                     
boundary value problem with λ being close to (Subscript[λ, 1]).  Let us use spectral Ritz-Galerkin method now. 

Corresponding potential:

(∫_0^1)(\(TraditionalForm\`-\+p\%1\)|u'((|)^p) - λ \(TraditionalForm\`-\+p\%1\)|u((|)^p) + ((f^T)+ a (Subscript[φ, 1]) ) u) dx ,

For p=2 k,  k ∈ ℕ we can handle this problem symbolically (at least partially).

We take

u(x) = 
            b1 sin(x) + b2 sin(2 x) + b3 sin(3 x) + 
              b4 sin(4 x)

as the finite dimensional approximation of u.




For p=4 and f(x)=sin(2 x), the \(integrand\ in\ the\ potential\ is\ computed\ \(\(by\)\(:\)\)\)


int = 
ReleaseHold[
Hold[(1/
                      p D[u[x], x]^p - λ/p (u[x])^p + 
                  f[x] u[x])] /. {u[x] -> 
                b1 Sin[x] + b2 Sin[2 x] + b3 Sin[3 x] + 
                  b4 Sin[4 x], 
f[x] -> 
                Sin[2 x]} 
]


\(1\/4\ \((b1\ Cos[x] + 2\ b2\ Cos[2\ x] + 3\ b3\ Cos[3\ x] + 4\ b4\ Cos[4\ x])\)\^4 + 
        Sin[2\ x]\ \((b1\ Sin[x] + b2\ Sin[2\ x] + b3\ Sin[3\ x] + 
              b4\ Sin[4\ x])\) - 
        1\/4\ λ\ \((b1\ Sin[x] + b2\ Sin[2\ x] + b3\ Sin[3\ x] + b4\ Sin[4\ x])\)\^4\)




Calculating the potential by direct integration:

Integrate[int,  {x,  0,  Pi}]  //  Timing


{27.139` Second,((3 (b1^4) π)/32) + ((b2 π)/2) + 
            (3/2) (b1^2) (b2^2) π + ((3 (b2^4) π)/2) + 
            (3/8) (b1^3) b3 π + (9/2) b1 (b2^2) b3 π + 
            (27/8) (b1^2) (b3^2) π + 
            (27/2) (b2^2) (b3^2) π + ((243 (b3^4) π)/32) + 
            3 (b1^2) b2 b4 π + 18 b1 b2 b3 b4 π + 
            27 b2 (b3^2) b4 π + 6 (b1^2) (b4^2) π + 
            24 (b2^2) (b4^2) π + 54 (b3^2) (b4^2) π + 
            24 (b4^4) π - (3/32) (b1^4) π λ - 
            (3/8) (b1^2) (b2^2) π λ - 
            (3/32) (b2^4) π λ + 
            (1/8) (b1^3) b3 π λ - 
            (3/8) b1 (b2^2) b3 π λ - 
            (3/8) (b1^2) (b3^2) π λ - 
            (3/8) (b2^2) (b3^2) π λ - 
            (3/32) (b3^4) π λ + 
            (3/8) (b1^2) b2 b4 π λ - 
            (3/4) b1 b2 b3 b4 π λ - 
            (3/8) b2 (b3^2) b4 π λ - 
            (3/8) (b1^2) (b4^2) π λ - 
            (3/8) (b2^2) (b4^2) π λ - 
            (3/8) (b3^2) (b4^2) π λ - 
            (3/32) (b4^4) π λ}



\(Tricky\ way\ of\ integration\ by\ using\ linearity\ of\ the\ \(\(integral\)\(:\)\)\)


Map[
Integrate[#,  {x,  0,  Pi}] &, 
          Expand[int]
]  //  Timing


{5.167` Second,((3 (b1^4) π)/32) + ((b2 π)/2) + 
            (3/2) (b1^2) (b2^2) π + ((3 (b2^4) π)/2) + 
            (3/8) (b1^3) b3 π + (9/2) b1 (b2^2) b3 π + 
            (27/8) (b1^2) (b3^2) π + 
            (27/2) (b2^2) (b3^2) π + ((243 (b3^4) π)/32) + 
            3 (b1^2) b2 b4 π + 18 b1 b2 b3 b4 π + 
            27 b2 (b3^2) b4 π + 6 (b1^2) (b4^2) π + 
            24 (b2^2) (b4^2) π + 54 (b3^2) (b4^2) π + 
            24 (b4^4) π - (3/32) (b1^4) π λ - 
            (3/8) (b1^2) (b2^2) π λ - 
            (3/32) (b2^4) π λ + 
            (1/8) (b1^3) b3 π λ - 
            (3/8) b1 (b2^2) b3 π λ - 
            (3/8) (b1^2) (b3^2) π λ - 
            (3/8) (b2^2) (b3^2) π λ - 
            (3/32) (b3^4) π λ + 
            (3/8) (b1^2) b2 b4 π λ - 
            (3/4) b1 b2 b3 b4 π λ - 
            (3/8) b2 (b3^2) b4 π λ - 
            (3/8) (b1^2) (b4^2) π λ - 
            (3/8) (b2^2) (b4^2) π λ - 
            (3/8) (b3^2) (b4^2) π λ - 
            (3/32) (b4^4) π λ}


Speed up:

27.139` /5.167`
5.252370814786143`

Here, we use the fact that expanded integrand has the head Plus :

Expand[int]

(1/4) (b1^4) Cos[x]^4 + 2 (b1^3) b2 Cos[x]^3 Cos[2 x] + 
        6 (b1^2) (b2^2) Cos[x]^2 Cos[2 x]^2 + 
        8 b1 (b2^3) Cos[x] Cos[2 x]^3 + 4 (b2^4) Cos[2 x]^4 + 
        3 (b1^3) b3 Cos[x]^3 Cos[3 x] + 
        18 (b1^2) b2 b3 Cos[x]^2 Cos[2 x] Cos[3 x] + 
        36 b1 (b2^2) b3 Cos[x] Cos[2 x]^2 Cos[3 x] + 
        24 (b2^3) b3 Cos[2 x]^3 Cos[3 x] + 
        (27/2) (b1^2) (b3^2) Cos[x]^2 Cos[3 x]^2 + 
        54 b1 b2 (b3^2) Cos[x] Cos[2 x] Cos[3 x]^2 + 
        54 (b2^2) (b3^2) Cos[2 x]^2 Cos[3 x]^2 + 
        27 b1 (b3^3) Cos[x] Cos[3 x]^3 + 
        54 b2 (b3^3) Cos[2 x] Cos[3 x]^3 + (81/4) (b3^4) Cos[3 x]^4 + 
        4 (b1^3) b4 Cos[x]^3 Cos[4 x] + 
        24 (b1^2) b2 b4 Cos[x]^2 Cos[2 x] Cos[4 x] + 
        48 b1 (b2^2) b4 Cos[x] Cos[2 x]^2 Cos[4 x] + 
        32 (b2^3) b4 Cos[2 x]^3 Cos[4 x] + 
        36 (b1^2) b3 b4 Cos[x]^2 Cos[3 x] Cos[4 x] + 
        144 b1 b2 b3 b4 Cos[x] Cos[2 x] Cos[3 x] Cos[4 x] + 
        144 (b2^2) b3 b4 Cos[2 x]^2 Cos[3 x] Cos[4 x] + 
        108 b1 (b3^2) b4 Cos[x] Cos[3 x]^2 Cos[4 x] + 
        216 b2 (b3^2) b4 Cos[2 x] Cos[3 x]^2 Cos[4 x] + 
        108 (b3^3) b4 Cos[3 x]^3 Cos[4 x] + 
        24 (b1^2) (b4^2) Cos[x]^2 Cos[4 x]^2 + 
        96 b1 b2 (b4^2) Cos[x] Cos[2 x] Cos[4 x]^2 + 
        96 (b2^2) (b4^2) Cos[2 x]^2 Cos[4 x]^2 + 
        144 b1 b3 (b4^2) Cos[x] Cos[3 x] Cos[4 x]^2 + 
        288 b2 b3 (b4^2) Cos[2 x] Cos[3 x] Cos[4 x]^2 + 
        216 (b3^2) (b4^2) Cos[3 x]^2 Cos[4 x]^2 + 
        64 b1 (b4^3) Cos[x] Cos[4 x]^3 + 
        128 b2 (b4^3) Cos[2 x] Cos[4 x]^3 + 
        192 b3 (b4^3) Cos[3 x] Cos[4 x]^3 + 64 (b4^4) Cos[4 x]^4 - 
        (1/4) (b1^4) λ Sin[x]^4 + b1 Sin[x] Sin[2 x] - 
        (b1^3) b2 λ Sin[x]^3 Sin[2 x] + b2 Sin[2 x]^2 - 
        (3/2) (b1^2) (b2^2) λ Sin[x]^2 Sin[2 x]^2 - 
        b1 (b2^3) λ Sin[x] Sin[2 x]^3 - 
        (1/4) (b2^4) λ Sin[2 x]^4 - 
        (b1^3) b3 λ Sin[x]^3 Sin[3 x] + 
        b3 Sin[2 x] Sin[3 x] - 
        3 (b1^2) b2 b3 λ Sin[x]^2 Sin[2 x] Sin[3 x] - 
        3 b1 (b2^2) b3 λ Sin[x] Sin[2 x]^2 Sin[3 x] - 
        (b2^3) b3 λ Sin[2 x]^3 Sin[3 x] - 
        (3/2) (b1^2) (b3^2) λ Sin[x]^2 Sin[3 x]^2 - 
        3 b1 b2 (b3^2) λ Sin[x] Sin[2 x] Sin[3 x]^2 - 
        (3/2) (b2^2) (b3^2) λ Sin[2 x]^2 Sin[3 x]^2 - 
        b1 (b3^3) λ Sin[x] Sin[3 x]^3 - 
        b2 (b3^3) λ Sin[2 x] Sin[3 x]^3 - 
        (1/4) (b3^4) λ Sin[3 x]^4 - 
        (b1^3) b4 λ Sin[x]^3 Sin[4 x] + 
        b4 Sin[2 x] Sin[4 x] - 
        3 (b1^2) b2 b4 λ Sin[x]^2 Sin[2 x] Sin[4 x] - 
        3 b1 (b2^2) b4 λ Sin[x] Sin[2 x]^2 Sin[4 x] - 
        (b2^3) b4 λ Sin[2 x]^3 Sin[4 x] - 
        3 (b1^2) b3 b4 λ Sin[x]^2 Sin[3 x] Sin[4 x] - 
        6 b1 b2 b3 b4 λ Sin[x] Sin[2 x] Sin[3 x] Sin[
            4 x] - 3 (b2^2) b3 b4 λ Sin[2 x]^2 Sin[3 x] Sin[
            4 x] - 3 b1 (b3^2) b4 λ Sin[x] Sin[3 x]^2 Sin[
            4 x] - 3 b2 (b3^2) b4 λ Sin[2 x] Sin[3 x]^2 Sin[
            4 x] - (b3^3) b4 λ Sin[3 x]^3 Sin[4 x] - 
        (3/2) (b1^2) (b4^2) λ Sin[x]^2 Sin[4 x]^2 - 
        3 b1 b2 (b4^2) λ Sin[x] Sin[2 x] Sin[4 x]^2 - 
        (3/2) (b2^2) (b4^2) λ Sin[2 x]^2 Sin[4 x]^2 - 
        3 b1 b3 (b4^2) λ Sin[x] Sin[3 x] Sin[4 x]^2 - 
        3 b2 b3 (b4^2) λ Sin[2 x] Sin[3 x] Sin[4 x]^2 - 
        (3/2) (b3^2) (b4^2) λ Sin[3 x]^2 Sin[4 x]^2 - 
        b1 (b4^3) λ Sin[x] Sin[4 x]^3 - 
        b2 (b4^3) λ Sin[2 x] Sin[4 x]^3 - 
        b3 (b4^3) λ Sin[3 x] Sin[4 x]^3 - 
        (1/4) (b4^4) λ Sin[4 x]^4


We can define:

pot[λ_,  p_,  f_, b1_, b2_, b3_, b4_,  b5_]:=Map[
Integrate[#,  {x, 0, Pi}] &,
\(Expand[
ReleaseHold[
Hold[\((1/p\ D[u[x], x]^p - λ/p\ \((u[x])\)^p + 
                          f[x]\ u[x])\)] /. 
                    u[x] -> 
                      b1\ Sin[x] + b2\ Sin[2\ x] + b3\ Sin[3\ x] + 
                        b4\ Sin[
                            4\ x]
]
]\)
]


Usage:

pot[λ,  3,  Sin[2  #] &,  b1, b2, b3, b4]



Parallelization is usefull here !


use ParallelMap instead of Map


To avoid a lot of communication, use partitioning.

Example of partitioning:

On a 5-node homogeneous cluster, do the following

ParallelMap[
Map[
(Integrate[#,  {x, 0, Pi}] &), #
]&,  
With[{numberOfNodes = 5},
   Partition[List@@\(Expand[
ReleaseHold[
Hold[\((1/p\ D[u[x], x]^p - λ/p\ \((u[x])\)^p + 
                                        f[x]\ u[x])\)] /. {u[x] -> 
                                      b1\ Sin[x] + b2\ Sin[2\ x] + 
                                        b3\ Sin[3\ x] + 
                                        b4\ Sin[4\ x], 
f[
                                        x] -> Sin[
                                        2\ x]}\ 
]
]\), numberOfNodes,numberOfNodes, {1, 1}, 0] // Transpose
]
]



More clear explanation to partitioning:

\(With[{numberOfNodes = 
            5}, 
Partition[
{int01, 
              int02, int03, int04, int05, 
int06, int07, 
              int08, int09, 
              int10, 
int11}, 
\ \ numberOfNodes, numberOfNodes, \ {1, 1}, \ 0]\ 
 // \ 
          Transpose
]\ \)
 // TableForm


({{processor 1 :      int01, int06, int11}, {processor 2 :      int02, int07, 0}, {processor 3 :      int03, int08, 0}, {processor 4  :       int04, int09, 0}, {processor 5  :       int05, int10, 0}})

\(\(\(
\)\(
\)\)\)
Equations to solve:

D[pot[λ,  3,  Sin[2  #] &,  b1, b2, b3, b4], b1] == 0
D[pot[λ,  3,  Sin[2  #] &,  b1, b2, b3, b4], b2] == 0
D[pot[λ,  3,  Sin[2  #] &,  b1, b2, b3, b4], b3] == 0
D[pot[λ,  3,  Sin[2  #] &,  b1, b2, b3, b4], b4] == 0


Using FindRoot (or Reduce if the number of variables is not very big) we can obtain bifurcation diagram in 3-D  ℝ×W where W is spanned by {Sin(x), Sin(2 x)}:

This method works also for partial differential equations (only prospective difficulty is to find basis functions in  (W^(1,p))(Ω)).

Parameter Idetification in Nonlinear Models

If a heat-exchanger tube array is subjected to a cross-flow  fluid-elastic instabilities occur for some range of velocities. Since  the resulting amplitudes of tube oscillations may become large causing damage of the heat exchanger,this phenomena of self-excitation has to be avoided whenever possible. Lacking satisfactory theoretical models  in aeroelasticity we are limited to semi-empirical models. On purpose, we employ a nonlinear model:

    m u'' + b u' + k u - (V^2) γ (u-y-c)(u-y)(u-y+c)=0
    y'+β y = u
    u(0) = (Subscript[u, 0]), u' (0) = (Subscript[v, 0]), y(0)=0

introduced by Tondl [Tondl, Quenching of Self-Excited Vibrations, Academia, Prague, 1991] in a slightly different context (here m, b, k are structural mass, damping, and stiffness, respectively, u is the deflection of the tube, y is an auxiliary variable, (V^2) the free-stream velocity,  c, β and γ are parameters, ).   

Realistic values of parameters c, β and γ are difficult to estimate theoretically. For that reason we have decided to estimate them using parameter identification from measured data. There were several problems we encountered in doing so.We mention the most significant ones:

(i) Ratio of the time needed to achieve either stable equilibrium or the limit cycle and the period of oscillations of the system is very high (about  30000). This forces to use  very small integration step when solving the system numerically (small in comparison with the time domain) and to use enormous amount of data to which the solution of the nonlinear problem has to be fitted.This leads to possible accumulation of rounding errors.

(ii) The objective function has many local critical points.To find the global minimum is very difficult.

(iii) Fluctuations of the velocity of the free stream during the measurement.

We have partially solved problem (i) by using higher precision arithmetic supported by Mathematica (if it would not have been enough we can use multishooting methods). This approach is extremely both memory and time consuming, however.To partially overcome problem (ii),we use global optimization algorithms based on differential evolution (genetic algorithm) to get near to the  prospective global minimum. Then we proceed by using standard local methods such as the Newton method to find the critical point. Finally,the problem (iii) was solved by considering the velocity of the free stream as a function of t.

Conclusion  -  improvements suggestions

In processing of our problems we often use shooting method, which is implemented using NDSolve accompanied  with FindRoot, ContourPlot or NMinimize in Mathematica. Since NDSolve is called with the same equation and different initial conditions many times within the shooting, we found NDSolve methods very useful (significant speed-up). In many of our problems parameters are present. However, NDSolve methods are unable to process new choices of parameters. For that reason we would like to ask developers to include this possibility in NDSolve methods together with initial conditions in one command if possible (one would reduce number of high-level calls resulting in increase of speed). Finaly, I would like to present an example where also options of NDSolve have to be changed in each step of computation of the bifurcation diagram.

- (|u'((|)^(p- 2))u')' - λ \(\(|\)\(u\)\( | \^\(p - 2\)\)\)u + Sin[u]= (f^T)+ a (Subscript[φ, 1]) ,               in  (0, 1);
                     u(0)=0,  u(1)=0 .

Analogously as above, we employ shooting method for obtaining bifurcation diagrams. It is known see e.g. [7,8] or references therein that  u/∥u∥ has its profile of the form Sin[π t] in [0,1] for ∥u∥-> +∞. For that reason,
we have to integrate initial value problem with highly oscillating term Sin[u]≈Sin[∥u∥ Sin[π t]]:

Plot[Sin[100*Sin[Pi t]],  {t, 0,  1}]

Plot[Sin[1000*Sin[Pi t]],  {t, 0,  1},  
      PlotPoints ->  3000]

In order to get reasonable results, we have to satisfy Nyquist criterion of the integrand sampling. Since the frequency of oscillations is strongly depending on ∥u∥ which can be estimated from u'(0), it would be nice to have possibility of  having NDSolve methods for changing MaxStepSize together with processing initial values, (and hopefully with parameters; see above). This would significantly speed-up our computations. 

Let me conclude this work advertising our project called CML (Central Mathematics Laboratory) where most of know-how presented in contribution will be accesible through webMathematica interface.

References to works of our research workgroup

[1] Benedikt, J., Girg, P., Takáč, P., The Fredholm Alternative for the p-Laplacian at higher eigenvalues in one dimension, in Preparation. 
[2] Čepička, J., "Numerical Experiments for Nonlinear Problems'',  Ph.D. Thesis, University of West Bohemia, Pilsen, 2001, pp. 97. 
[3] Čepička, J.,  Drábek, P., Girg, P., Quasilinear  Boundary Value Problems: Existence and Multiplicity Results, in 
"Variational Methods: Open Problems, Recent Progress, and Numerical Algorithms", John M. Neuberger, Northern Arizona University, Editor, Contemporary Mathematics Series, AMS, 2004, 285 pp., Softcover, ISBN 0-8218-3339-1.
[4] Čepička, J., Drábek, P., Girg, P., Open problems related to the p-Laplacian operator,   
Bol. Soc. Esp. Mat. Apl. 29(2004), pp. 13-34. 
[5] Del Toro, P., Girg, P.,  Roca, F., Dirichlet Problem with Nonlinearity Depending only on the Derivative,  Appl. Math. Lett. 16(2003), pp. 7-12. 
[6] Drábek, P., Girg, P.,  Roca, F., On the Range Properties of Certain Landesman-Lazer-type Problem, J. Math. Anal. Appl.  257(2001), pp. 131-140.  
[7] Drábek, P.,  Girg, P., Takáč, P., Ulm, M., The Fredholm alternative for the p-Laplacian: bifurcation
from infinity, existence and multiplicity of solutions, Indiana Univ. Math. J., 53 (2) (2004), pp. 433-482.
[8] Drábek, P., Girg, P., Takáč, P., Bounded  perturbations of homogeneous quasilinear operators using bifurcations from infinity, to appear in J. Differential Equations,  accepted 2003, electronic version available on-line via www.elsevier.com (2004).
[9] Frisch J., Personal Communication (permanent).
[10] Girg, P., Some Notes on Aeroelastics of Tube Arrays, Technical Report, University of West Bohemia, Pilsen 2004. 
[11] Girg, P., Roca, F., Villegas, S., Semilinear Sturm Liouville Problem with Periodic Nonlinearity, to appear in Nonlinear Analysis T.M.A., accepted 2000.
[12] Girg, P., Takáč, P., Bifurcations of Positive and Negative Continua in Quasilinear Elliptic Eigenvalue Problems,
in Preparation.
[13] Girg, P., Mathematical Model of Heat-Exchanger Tube and Parameter Identification, in "Interaction of Dynamic Systems with Surroundings and Systems with Feedbacsk", I. Zolotarev, Institute of Thermomechanics AS CR, Editor, Prague, pp. 17-22, Softcover, ISBN  80-85918-91-9.

References
[1]  Mathematica Book 3-5
[2]  NDSolve Advanced Documentation
[3]  Parallel Computing Toolkit Documentation
[4]  All references cited in works in References to works of our research  workgroup

www . cam . zcu . cz

drift . cam . zcu . cz

Acknowledgement

All expanses related to my participaton at Wolfram Developer Conference 2004 and preparation of this contribution have been covered by the grant 1N04078 of the Ministry of Education, Youth and Sports of the Czech Republic (sponsor of CML). Since numerical experiments are important tool in my theoretical research in various topics in nonlinear analysis,  research presented in this contribution has been supported from several sources depending on subject
(in chronological ordering):
by the Grant number 201/03/0671 (quasilinear differential equations; nonlinear stationary diffusion); 
by the Grant number 101/02/1225 (nonlinear models of heat exchanger tubes);
by the Humboldt Foundation during my 12-month fellowship 2003/04 (quasilinear differential equations with periodic nonlinearity).

 Finally, I would like to thank to Wolfram Research staff for organizing very nice conference and especially for possibility of consulting and solving our problems. My special thanks go to Joy Costa, Rob Knapp ,  Roman Maeder, and Michael Trott .