Numerical Analysis

Archive for December, 2009

matrix-vector multiplication in LINPACK and LAPACK

I would like to know what storage schemes are used to store
matrices in LINPACK and LAPACK.

Also, how to do matrix-vector multiplication in those storage schemes.

Thanks in advance.

–Ben Lam
l…@research.cs.orst.edu

.
posted by admin in Uncategorized and have No Comments

Automatic Differentiation : Help needed

G’day to all num-analysis netters,

I’ve just been recommended to seek assistance from you guys regarding
the subject of automatic differentiation to solve a problem regarding
symbolic programming in Fortran.

First of all, here’s the background of the problem that I have :

I have at present a Fortran program which is numerically intensive
at the moment. I need to add a section to that program which is
unfortunately symbolically intensive. The requirements of the
symbolically intensive portion is to find the partial derivatives
of a multi-variable scalar function s(x,y,z) which is given by

    s(x,y,z) =  w(T) . Q . w  

 where

    (T) denotes transpose of a vector or matrix

     w is an input vector whose elements consists of constants and
     variables –
     eg : w = [2 , 2x+y, x^2, y^4, 3x , 4x-z, 4, 5.6 ](T)

     Q is an input matrix whose elements are similar to the elements
     in w, ie , they can be symbolic elements (ie in terms of variables
     x, y or z) or purely numerical entries. The symbolic elements in
     Q are limited to integer powers in x, y or z .

 Question is , how do you compute

   d s(x,y,z)          d s(x,y,z)           d s(x,y,z)
   ———-     ,    ———-     and   ———-       ?
      d x                 d y                  d z

I guess to do the above in mathematica or maple is simple enough :

a) Define w and Q and multiply out w(T).Q.w which basically makes s(x,y,z)
   to be a polynomial function in x,y and z
b) Then collect the terms and apply the partial differentiation function
   on s(x,y,z)

Well, the unfortunate part is that all these operations have to be done
within Fortran 77.

Someone have recommended me to use ‘automatic differentiation’ . I haven’t
got the faintest idea what that is and how could it be applied to solve
the above problem in Fortran. I am sure that someone out there has seen
the above problem before and hopefully have found a way out of this mess
using automatic differentiation.

If there is anyone out there who knows what is automatic differentiation
or have come across references on this subject, I would really appreciate
if you could tell me about it and even better, show me how I can apply
it to solve the above symbolic problem in Fortran.

A thousand thanx to advance to anyone willing to help me out !

I can be reached at :

——————————————————————————
                 Ian Thng  - Forever Trying to be Young  
                 Email:t…@atri.curtin.edu.au
                 Adaptive Signal Processing Laboratory
                 Australian Telecommunications Research Institute
                 Curtin University of Technology
                 Ph :  +61-9-3513270
                 Fax : +61-9-3513244
                 Bothered with seat belts ? Ride a Harley Mate !
——————————————————————————

posted by admin in Uncategorized and have Comments (4)

Questions re convolution in Num. Rec. in C

Hello,
        I have been trying to use the code from the Numerical Recipes
in C book.  I’m not sure which edition (there are two aren’t there?) I
copied the code from.  However, I am not getting correct values for the
convolution.  Has anyone had any experience with these routines
(convolution and the fft stuff)?  Are there any bugs in the code?  I’ve
checked mine with the book a number of times and feel 95% confident that
my code is consistent with the book.  
        Any help is appreciated.  This convolution is a small part in a
bigger program for my research and I’ve been bogged on it for many weeks
now.  Someone would be much in my debt if they can shed light here.

charles d. kincaid

p.s.  btw, am I alone in believing that the documentation for the code
is very poor.  They do a good job explaining the theory behind the code
(for fft stuff at least), but then they use bad variables and don’t
explain why the code works.  It could be just me.  

chuck

charles d. kincaid              STATS:  (904) 392-1941
stat…@stat.ufl.edu          CIRCA:  (904) 392-2007
‘Damn fine coffee…and hot, too!’

posted by admin in Uncategorized and have No Comments

Finite elements and Navier-Stokes equations

Hi, I working on devloping a finite-element model to model the wind
driven barotropic oceanic circulation. My background, as well of that
of my supervisor) is mainly oceanography, and I’ve reached the stage
where I have several questions where somebody with more familiarity
with the field may be able to point me towards some useful previous
work and references.
   My basic equation of interest is

    epsilon*(grad**4)(psi)+d(psi)/dx+J(psi,grad**2(psi))=forcing

where: epsilon is a small parameter, 6.0e-6
psi is the barotropic streamfunction that I wish to solve for
J(psi,grad**2(psi)) is the Jacobian
d/dx is the partial derivative (just can’t indicate it elseways)

   Now, as far as I can see, there are two ways to solve this:
1) Direct solution of the fourth order equation. Even after
integrating by parts I will still have to take second derivatives
of my basis functions, so I need to use C1 continuity elements.
I wish to use triangular elements. I’ve read about an 18 node
incomplete quintic element that has been used (e.g. Olson & Tuann),
but I have to admit I can’t get this to work. I was wondering if
anyone had some experience with this (or some detailed papers)
about implementing this element.
2) Solving using the coupled Streamfunction-Vorticity Equations.
Here, the vorticity is defined as the Laplacian of the streamfunction
giving my a system of two coupled p.d.e.’s. Using Green’s identity
I can reduce the order of each equation to first, so i can use
simple C0 continuity elements. However, the problem with this comes
in with the boundary conditions. My boundary conditions are
psi=0 on the boundary and d(psi)/dn=0 on the boundary (no slip).
So i have two boundary conditions on the streamfunction and none on
vorticity, which is the problem. If I want to use slip conditions
(vorticity=0 on boundary), I can get this to work with no problem.
   But getting a boundary condition on vorticity with no slip is my
problem. I’ve read about several ways to do this.
a – Expand grad**2(psi)=w {w is my vorticity} in a finite difference
expansion and get w on the boundary in terms of psi normal to it,
but this would require me to have a node normal to all boundary nodes,
and i think this method isn’t as numerically accurate either.
b – Let the vorticity on the boundary be represented by some constant
lambda, where lambda has to be found previously. This is based on some
work by Glowinski and Pironneau, I beleive, but again I haven’t been
able to implement it in practice.
c – A seemingly direct method that directly solves the equations, as
outlined in papers by Campion-Renson and Crochet, Dhatt et al., and
Mizukami (among others), all found in the International Journal of
Numerical Methods in Engineering. However, when I try to implement
this, I get what seem like instabilities in my solution, with alternating
positive and negative gyres. I was wondering if anyone had seen this
before and maybe new a possible fix (increased resolution, different
basis functions, ???).

   Anyways, I have probably explained this in too much detail but have
spent a lot of time on this, and don’t seem to be getting anywhere on
this last little bit, so I was wondering if anyone had any experience
with problems such as this or new some more good references.
   Also, I would be interested in hearing which technique is preferred,
solving the direct fourth order biharmonic equation or using the
streamfunction-vorticity technique, and why.
   Thanks

Paul Myers
School of Earth and Ocean Sciences,
P.O. Box 1700, University of Victoria,
Victoria, B.C., Canada
(604) 721-6193
my…@ocean.seaoar.uvic.ca

posted by admin in Uncategorized and have No Comments

Re: Bessel Fcns of 2nd kind of integer order

ba…@ramsey.cs.laurentian.ca (Barry G. Adams) writes:

>The IMSL special function library doesn’t seem to have a routine
>for calculating
>                       Y  (z)
>                         n
>where n is an integer and z = x + iy is a complex number. Does
>anyone know how to do this or do you know of an available routine?

Algorithm 644 from Transactions on Math Software is a massive collection
of portable Bessel function codes by Donald Amos, and includes such      
routines. You can get the whole thing by sending email to

     net…@ornl.gov

with the message

    send big644 from toms
or  send 644r from toms

The response will be LARGE (over 600K) and will arrive in pieces.

To be more selective about which routines you want from the set, send
the message

    send index for amos

You’ll receive a list of items available in the amos collection; you can
then order individual items with messages like

    send x from amos

Stan Kerr    
Computing & Communications Services Office, U of Illinois/Urbana
Phone: 217-333-5217  Email: stank…@uiuc.edu  

posted by admin in Uncategorized and have Comments (2)

Conference on Graduate Programs in the Applied Mathematical Sciences

             Conference on Graduate Programs in the Applied
                       Mathematical Sciences II
                         April  16-18, 1993
                         Clemson University
                            Clemson, SC

Current challenges for graduate education in classical analysis,
computational and discrete mathematics, operations research and
statistics will be discussed at this conference hosted by the
Department of Mathematical Sciences at Clemson University.

Specific conference objectives are to review the current status of
graduate education in the applied mathematical sciences; to review
the role of recent developments in the mathematical sciences in
the solution of modern technological problems and the challenges
and opportunities these developments present to graduate
education; to obtain perspectives on current industrial/government
employment opportunities and placement; and to publicize
innovative graduate programs to the professional community.

Participants will be selected to provide broad representation from
all areas of the mathematical sciences, including individuals from
academic, governmental and industrial institutions. The program
will consist of a number of panel, audience, and group
discussions.  This format necessitates a limited number of
participants.  A conference proceedings will be published.

Preliminary Program
Friday April 16
Morning
   I. Graduate Mathematical Sciences Education: An Historical Perspective.
  II. Current Status of Applied Mathematical Sciences Graduate Education
      in the U.S.
 III. Challenges for Graduate Education in the Mathematical Sciences.

Afternoon
  IV. Curricula for the 90′s and Beyond.
   V. Curriculum Innovations.

Saturday April 17
Morning
  VI. Computation in the Curriculum.  
 VII. Employment Opportunities for the 90′s.  

Afternoon
VIII. Non-Academic Careers for MS and Ph.D. Students.
  XI. Group Discussions.
      1. Recruiting and Advising Graduate Students.
      2. Department Administration, University Support, Rewards.
      3. Undergraduate Preparation.

Sunday April 18
  X.  Summary and Closing Remarks.

Invited participants include: John Burns, Don Bushaw, Ron Douglas, James
Dyer, Bob Easterling, Avner Friedman, Carl Harris, Robert Hogg, Bill Lucas,
Robert Lundegard, Tom Magnanti, Joyce McLaughlin, Anna Nagurney,
Fred Roberts, Mary Wheeler, Doug Zahn.    

This conference is supported by a grant from the National Security
Agency and additional funding from NSF is expected.  For further
information contact, R. E. Fennell, Department of Mathematical
Sciences, Clemson University, Clemson SC 29634-1907, (803) 656-3257,
email c…@math.clemson.edu  

                Matthew Saltzman
                Clemson University Math Sciences
                m…@clemson.edu

posted by admin in Uncategorized and have No Comments

Problem: Computing path length on cubic spline?

Hi Netters,

  I’ve stumbled on what appears to be a surprisingly intractable
problem.  Let us say that you are using a set of cubic spline
functions to define a path through a set of points.  Now you
want to know the total length of this path.  Simple, right?
For each spline function,

         y = A + Bx + Cx^2 + Dx^3
so
        dy = [b + 2cx + + 3Dx^2]

so the path length is described by a differential equation of the form

        dl = SQRoot[ dx^2 + dy^2] = SQRoot[Q + Rx + Sx^2 + Tx^3 + Ux^4]dx

  I have come to despair of finding an indefinite or definite
closed-form solution to the integral of dl.  Is there no
alternative to integrating these spline functions numerically? (Ouch!)
I thought that cubic splines were pretty straightforward, but
suddenly I’m in over my head.  If anyone has encountered this
problem before, or cares to make a suggestion, I would be grateful.

Sean Burke

Sean Burke                      s…@netcom.com

posted by admin in Uncategorized and have Comment (1)

Software Package(s) for PDEs

Hello Netters,

I am looking for a software package (preferably in C)
that allows for the numerical solution of the Euler
equations (i.e., the Navier-Stokes equations with the
viscosity terms = 0).  The system is in 3-D (i.e.,
x,y,z coordinates plus time, t). I am contemplating
running this package either on PC or on a SUN Workstation.
Therefore, if executables are not available, I would
adapt source code to them.

Please respond to k…@hogpa.ho.att.com.

Thank you very much in advance,

John C. Kiss
k…@hogpa.ho.att.com

posted by admin in Uncategorized and have No Comments

Solving Navier-Stokes Eqn using Penalty Function Method

Hi,
I am solving Navier-Stokes Eqn. using FEM (Penalty-function method) right
now. I have read lots of papers, but nobody mentioned how they handeled
the force vector. The force vector is

        F1=line or suface integration of(tx)
        F2=line or suface integration of(ty)
where
        tx= 2.0*mu*(du/dx)*nx+mu*((du/dy)+(dv/dx))*ny-p*nx
        ty= mu*((du/dy)+(dv/dx))*nx+2.0*mu*(du/dx)*nx-p*nx
and
        p=-(penalty number)*(du/dx+dv/dy)

I am now using iteration method to handel this, ie first I assume a velocity
profile, secondly calculate the force vector and then sove for the velocity.
After this , I use the vel. got for the next iteration. But it seems that
the result does not converge.

Is there anyone can give me some help?

Regards.

lixin

email: lz…@uoft02.utoledo.edu      

posted by admin in Uncategorized and have Comment (1)

<None>

Hi,
I am solving Navier-Stokes Eqn. using FEM (Penalty-function method) right
now. I have read lots of papers, but nobody mentioned how they handeled
the force vector. The force vector is

        F1=line or suface integration of(tx)
        F2=line or suface integration of(ty)
where
        tx= 2.0*mu*(du/dx)*nx+mu*((du/dy)+(dv/dx))*ny-p*nx
        ty= mu*((du/dy)+(dv/dx))*nx+2.0*mu*(du/dx)*nx-p*nx
and
        p=-(penalty number)*(du/dx+dv/dy)

I am now using iteration method to handel this, ie first I assume a velocity
profile, secondly calculate the force vector and then sove for the velocity.
After this , I use the vel. got for the next iteration. But it seems that
the result does not converge.

Is there anyone can give me some help?

Regards.

lixin

email: lz…@uoft02.utoledo.edu      

posted by admin in Uncategorized and have No Comments