Numerical Analysis

Archive for August, 2010

discrete nonlinear Burgers Eq.

Can someone help me understand how to set up to discretize the
velocity term in non-linear Burgers equation.  The equation is Ut/dt +
Uxx/2∙dx = nu∙Uxx/dxx.  Applyed general form dU/dt + dF/dx = 0 where
F=U^2/2 – nu∙dU/dx. Literature has "change order of differentiation"
of "F" in this form to get the flux Jacobian "A" and solve for the
Taylor series?   I have an example of linear case, but non-linear
terms are getting confusing.


posted by admin in Uncategorized and have Comments (3)

condition number issue


I would like to show that for a square matrix A, the following relation
holds approximately (at least for large values of k(?)):

cond(A^(k+1)) / cond(A^k) \approx EWmax/EWmin

with cond(.) the standard condition number using the 2-norm (hence the
largest divided by the smallest singular value of A) and EWmax the
largest eigenvalue of A (in magnitude) and EWmin the smallest eigenvalue
of A in magnitude.

It is probably trivial, but how do I prove this in an easy manner? Are
there any conditions for which the above equation holds more tightly, or
where it fails? Is it crucial that k is large?

Thanks in advance,

posted by admin in Uncategorized and have Comments (7)

What about roots ?

If k, n  are in {1,2,…} and (z)_k: =z(z+1)…(z+k-1) ,

what about the roots  of polynomial equation

(x)_n +(x+n+1)_n = 0   ?

posted by admin in Uncategorized and have Comments (2)

complex reciprocal Gamma around negative integers

I am looking for code (in C(++) ) to evaluate 1/Gamma(z), z = -k + eps,
eps complex and small, say |esp| <= 1E-4 and smaller with k a natural
number (not too large, say <=8(?) ).

Currently I tried recursion (with product = exp(log(sum)) ) to reduce
to 1/Gamma(1 + eps), but can not achieve absolute accuracy  ~ 1E-14, I
loose too much (using Maple for a cross check) and get ~ 1E-11 absolute
and even worse for relative exactness (in not too idiotic cases …).

Is there some *sound* code available or do I expect too much for double
precision? If it covers the right half plane as well it would be really
great (ok, Xmas is gone, but soon we have Eastern).

posted by admin in Uncategorized and have Comments (17)

Multiscale Methods and Statistics: A Productive Marriage

Call for submissions on

Following the successful 5th Graybill Conference on "Multiscale
Methods and Statistics: A Productive Marriage" organized by Department
of Statistics, Colorado State University in June 2006, we are inviting
submissions on the same theme for a special issue of Statistica
Sinica.  Submissions are NOT limited to papers presented at the
conference.  The goal of publishing this theme topic is to highlight
the recent developments and advantages of employing both multiscale
and statistical methods for solving various methodological, scientific
and engineering problems.

Submissions preferably will provide either methodological or
theoretical advances in the interface of multiscale methods and
statistics, or demonstrate applications of such techniques in sciences
and technologies. Priority will be given to papers that are assessed
to be, statistically and scientifically speaking, most innovative,
comprehensive, and reaching a wider readership.

Papers submitted to the special issue will be reviewed according to
regular procedures similar to other submissions to the journal.
Accepted papers will appear in a single issue of Statistica Sinica,
scheduled for 2009. Please contact Professor Yazhen Wang or any of the
special issue editors listed below for questions on the suitability of
your paper(s). Submissions must be made online through the journal
site at


Please use the LaTeX article template, also available at the above
site, for preparing your manuscript submission. The scheduled DEADLINE
for submission is December 31, 2007. Authors wishing to receive email
alerts as the deadline approaches may contact the editorial assistant,
Karen Li (email: karen at

Special Issue Editors:

T. Tony Cai, U. of Pennsylvania (tcai AT
Thomas Lee, The Chinese U. of Hong Kong (tlee AT
Yazhen Wang, U. of Connecticut (yzwang AT
Patrick J. Wolfe, Harvard U (patrick AT

Related Workshop


posted by admin in Uncategorized and have No Comments

Meshless to solve second order differential system of equations

Dear all,

Using Maple I am able to compute a general solution to my system of
ODE; does anyone know an easy way to find the particular solution to
the same problem using meshless methods, i.e. obtaining a solution
that verifies prescribed boundary conditions.

Here’s the Maple worksheet with my differential equations system


- Hide quoted text — Show quoted text -

> with(LinearAlgebra):
> with(DEtools):workdir:="D:/Science/PhD@work/3Dmodel/kirchhoff/general/BOX8LL":
> currentdir(workdir);
> read "modelRigid_kirchhoff.m":
> read "modelRigid_modesB2.m":
> Yc:=-B/2:
> Zc:=-H/2:
> alias(alpha=B^4+6*H*B^3-5*H^2*B^2+6*B*H^3+H^4);
> alias(beta=9*E*B^4*G*H^4+18*E*B^3*G*H^3*t^2+3*E*B^2*G*H^4*t^2+3*E*B^4*G*H^2*t^2+8*E*B^2*G*H^2*t^4+4*E*B*G*H^3*t^4+4*E*B^3*G*H*t^4);

> T_modesGlobal_t:=transpose(T_modesGlobal):
> K0g2d:=map(simplify,evalm(T_modesGlobal_t &* K0g2 &* T_modesGlobal));
> subs(Yc=-B/2,Zc=-H/2,evalm(K0g2d));
> nDim:=rowdim(K0g2d);
> K1g2d:=map(simplify,evalm(T_modesGlobal_t &* K1g2 &* T_modesGlobal));
> subs(Yc=-B/2,Zc=-H/2,evalm(K1g2d));

> K2g2d:=map(simplify,evalm(T_modesGlobal_t &* K2g2 &* T_modesGlobal));
> subs(Yc=-B/2,Zc=-H/2,evalm(K1g2d));
> d0:=vector(11,[q1(x),q2(x),q3(x),q4(x),q5(x),q6(x),q7(x),q8(x),V(x),W(x),Theta(x)]);
> d1:=map(diff,d0,x);
> d2:=map(diff,d0,x,x);
> EQD:=map(simplify,evalm(K2g2d &* d2 + K1g2d &* d1 + K0g2d &* d0)):
> EQD_num:=map(evalf,subs(H=2,B=1,t=0.01,E=200,G=87,evalm(EQD)));
> #for i from 1 to nDim do eq[i]:=EQD[i]=0: end do:
> for i from 1 to nDim do eq[i]:=EQD_num[i]=0: end do:
> sys_eqd:=[seq(eq[i],i=1..nDim)]:
> sol:=map(evalf,dsolve(sys_eqd));
> subs(sol,q1(x));
> subs(sol,q2(x));
> subs(sol,q3(x));
> Q0g2d:=evalm(T_modesGlobal_t &* Q0g2 &* T_modesGlobal);
> Q1g2d:=map(simplify,evalm(T_modesGlobal_t &* Q1g2 &* T_modesGlobal));

> dd1:=map(convert,d1,D);
> dd0:=map(convert,d0,D);
> Boundary_Terms:=evalm(Q1g2d &* dd1 + Q0g2d &* dd0);
> P:=vector(11,[N1,N2,N3,N4,N5,N6,N7,N8,Vy,Vz,T]);
> Pd:=map(simplify,evalm(T_modesGlobal_t &* P));
> for i from 1 to nDim do eq_boundary_static[i]:=Boundary_Terms[i]-Pd[i]=0; end do:
> EQBStatic:=[seq(eq_boundary_static[i],i=1..nDim)]:
> EQBStatic_num:=subs(x=0,H=2,B=1,t=0.01,E=200,G=87,N5=1000,N1=0,N2=0,N3=0,N4=0,N5=0,N6=0,N7=0,N8=0,Vy=0,Vz=0,T=0,EQBStatic);
> d_prescribed:=vector(11,[U1,U2,U3,U4,U5,U6,U7,U8,V_p,W_p,Theta_p]);
> EQBKinematic:=[seq(d0[i]-d_prescribed[i]=0,i=1..nDim)];

> EQBKinematic_num:=subs(x=10,H=2,B=1,t=0.01,E=200,G=87,U1=0,U2=0,U3=0,U4=0,U5=0,U6=0,U7=0,U8=0,V_p=0,W_p=0,Theta_p=0,EQBKinematic);
> sol_num:=dsolve({op(sys_eqd),op(EQBStatic_num),op(EQBKinematic_num)},{q1(x),q2(x),q3(x),q4(x),q5(x),q6(x),q7(x),q8(x),V(x),W(x),Theta(x)},numeric);

Is there an easy way to implement a meshless method to solve this out?

posted by admin in Uncategorized and have No Comments

JOB: Mathematician/contractor

Job Description

Title: Mathematician (Consultant/Contractor)

GFT Group seeks a highly creative individual with solid skills and
experience in advanced mathematics for a speech recognition project.
The project will develop new concepts and advanced mathematics that
substantially improve the accuracy of current speech recognition
technology.  This is a challenging project that requires unusual
skills.  It involves a close integration of advanced theory and
practice not found in most commercial algorithm projects and academic
research projects.  The project offers the opportunity to perform
fundamental scientific research with substantial immediate practical
benefits and applications.

GFT Group is developing a new speech recognition engine with superior
accuracy compared to current speech recognition engines such as Dragon
Naturally Speaking, ViaVoice, Microsoft Speech, and the open-source
SPHINX speech recognition engine.  The engine should enable or enhance
rapid jumping to targets such as files, web sites, PowerPoint slides,
submenu selections and so forth on computers by simple voice command,
dictation of documents, and hands-free operation of computers,
cell-phones, automobile accessories and many other devices.  Our
is to achieve 100% accurate speaker-independent phoneme recognition in
the presence of typical background sounds such as car noises that do
not impair human speech recognition.  We do not expect to solve the
homonym resolution problem for general unstructured human speech. The
effective speech recognition accuracy, the word error rate, of the
engine will be determined by the frequency of homonyms and
near-homonyms, words and phrases that sometimes sound the same, in the
recognized speech.

The Mathematician will help translate advanced concepts in human
to specific mathematical formulas that can be tested on human speech
data and, if successful, converted quickly to software for a real-time
commercial speech recognition engine written in a portable compiled
language such as ANSI C.  The engine will include a Microsoft Speech
compatible wrapper.

The task is similar to the inference of mathematical formulas such as
differential equations from experimental data and from concepts
expressed in words, pictures, and rough mathematical formulas.  It may
resemble, for example, the translation of Michael Faraday’s ideas
about electricity and magnetism from the words and pictures that
Faraday used to a set of new differential equations by James Clerk
Maxwell.  Experience with this process is most valuable for this
position but is not a requirement.

A strong knowledge of the physical processes and/or visual
representations corresponding to individual terms and factors within
terms in differential equations and other mathematical formulas should
be helpful. A strong knowledge of English verbal descriptions, words,
and phrases used for these physical processes and/or visual
representations should also be helpful.  This should make it easier to
identify the mathematics — for example construct a new differential
equation — corresponding to features of data and advanced concepts
expressed in words and pictures.

Experience in the following areas may be helpful, but is in no way a
specific requirement of the position:

1. Non-linear differential equations
2. Classical invariant theory
3. Differential geometry and tensor analysis
4. Hidden Markov Model speech recognition methods
5. Adaptive signal processing
6. Music theory and practice
7. Acoustics of Speech Communication
8. Implementation of advanced mathematical algorithms in C or similar
9. Mechanical modelling, finite element analysis, and simulation of
elastic materials.

The Mathematician must be able to think creatively, "outside the
box", and have a solid foundation in advanced mathematics including
the ability to learn new areas rapidly as needed.  An advanced degree
in mathematics, applied mathematics, or theoretical physics may be
helpful, but is not required.  This is primarily a pencil and paper
task; computer skills (for example, knowledge of a computer algebra
system such as Mathematica) may be helpful but are not required.

The Mathematician will be an independent contractor and will need to
sign a non-disclosure agreement to perform the task.  The
will need to work closely with a scientist from GFT Group in a
collegial style. This is not implementing someone else’s ideas but
requires shared creativity by skilled experts in complementary fields.
The opportunity to share license fees from the technology and inventor
status (on patent applications for example) may exist.

GFT Group is a private research and development and contracting
organization with several products in the speech recognition field.
For more information, see

Please submit a cover letter and a résumé and/or curriculum vitae in
the body of a single plain text e-mail message.  A statement of
research interests and research philosophy may be helpful, but is not
required.  Please send plain text only.  No attachments, hypertext, or
active content.  Please include the word "Mathematician" in e-mail
subject line.  Please send to:

John F. McGowan, Ph.D.
President, Research and Development Division
GFT Group Inc.

posted by admin in Uncategorized and have No Comments

Symposium "Computational Methods in Image Analysis" within the USNCCM IX Congress – Announce & Call for Papers


(Apologies for cross-posting)

Symposium "Computational Methods in Image Analysis"
National Congress on Computational Mechanics (USNCCM IX)
San Francisco, CA, USA, July 22 – 26 2007

We would appreciate if you could distribute this information by your
colleagues and co-workers.


Dear Colleague,

Within the National Congress on Computational Mechanics (USNCCM IX),
to be held in San Francisco, CA, USA, July 22 – 26 2007, we are
organizing the Symposium "Computational Methods in Image Analysis".
Examples of some topics that will be considered in that symposium are:
Image acquisition, Image processing and analysis, Image segmentation,
3D Vision, Motion analysis, Pattern recognition, Objects recognition,
Medical imaging and Tools and applications.
Due to your research activities in those fields, we would like to
invite you to submit an abstract and participate in the symposium
"Computational Methods in Image Analysis".

For instructions and submission, please access to the congress website
Please note, when submitting your abstract you should select
"Computational Methods in Image Analysis".

Important dates:

- Deadline for abstract submissions: May 1, 2007
– Final selection of abstracts: May 15, 2007
– Deadline for early registration: May 22, 2007
– Deadline for CD-ready abstracts: June 1, 2007
– Deadline for regular registration: July 15, 2007

We would appreciate if you could distribute this information by your
colleagues and co-workers.

Kind regards,

João Manuel R. S. Tavares
Renato Natal Jorge
Yongjie Zhang
Dinggang Shen
(Symposium organizers)

posted by admin in Uncategorized and have No Comments

a problem in numerical matrix inversion


     I am using some external package to deal with matrix operations,
mostly matrix inversion and trasposition. The dimentions of the
matrices used are about 29×19,19×19 and 29×29.

     During the calculation, I noticed an apparent error of
inverion of a 19×19 matrix. Denote this matrix as KK, U=KK^ -1, I
found the product of U and KK is not equivalent to unit matrix! This
apparently violate the definition of inversion. After some
I found the origin of the error is that the ratio of max/min
for this matrix is 7.4368432669e+016,Maybe this exceed the of
of the machine.

     Is there any tricks for me to be able to deal with this matrix

posted by admin in Uncategorized and have Comments (3)

Problem with DGESVD – Singular Value Decomposition (LAPACK)

For the matrix


calling DGESVD with argument JOBU=’N’ and JOBVT=’N’ the Compaq Visual
Fortran compiler goes into an infinite loop, while the g95 compilers
returns division by zero.

calling DGESVD with arguments JOBU=’A’ and JOBVT=’N’ the Compaq Visual
Fortran compiler returns log: -SING error, while the g95 compiler
returns the correct singular values that also also match MATLAB
function svd results.

Anyone has experience with the DGESVD routine or any ideas?


posted by admin in Uncategorized and have No Comments