Items tagged with numerical


I'm back from presenting work in the "23rd Conference on Applications of Computer Algebra -2017" . It was a very interesting event. This fifth presentation, about "The Appell doubly hypergeometric functions", describes a very recent project I've been working at Maple, i.e. the very first complete computational implementation of the Appell doubly hypergeometric functions. This work appeared in Maple 2017. These functions have a tremendous potential in that, at the same time, they have a myriad of properties, and include as particular cases most of the existing mathematical language, and so they have obvious applications in integration, differential equations, and applied mathematics all around. I think these will be the functions of this XXI century, analogously to what happened with hypergeometric functions in the previous century.

At the end, there is a link to the presentation worksheet, with which one could open the sections and reproduce the presentation examples.

The four double-hypergeometric Appell functions,

a complete implementation in a computer algebra system


Edgardo S. Cheb-Terrab

Physics, Differential Equations and Mathematical Functions, Maplesoft


The four multi-parameter Appell functions, AppellF1 , AppellF2 , AppellF3  and AppellF4  are doubly hypergeometric functions that include as particular cases the 2F1 hypergeometric  and some cases of the MeijerG  function, and with them most of the known functions of mathematical physics. Appell functions have been popping up with increasing frequency in applications in quantum mechanics, molecular physics, and general relativity. In this talk, a full implementation of these functions in the Maple computer algebra system, including, for the first time, their numerical evaluation over the whole complex plane, is presented, with details about the symbolic and numerical strategies used.

Appell Functions (symbolic)



The main references:


P. Appel, J.Kamke de Feriet, "Fonctions hypergeometriques et Hyperspheriques", 1926


H. Srivastava, P.W. Karlsson, "Multiple Gaussian Hypergeometric Series", 1985


24 papers in the literature, ranging from 1882 to 2015


Definition and Symmetries


Polynomial and Singular Cases


Single Power Series with Hypergeometric Coefficients


Analytic Extension from the Appell Series to the Appell Functions


Euler-Type and Contiguity Identities


Appell Differential Equations


Putting all together


Problem: some formulas in the literature are wrong or miss the conditions indicating when are they valid (exchange with the Mathematics director of the DLMF - NIST)


Appell Functions (numeric)






Compute these Appell functions over the whole complex plane


Considering that this is a research problem, implement different methods and flexible optional arguments to allow for:

a) comparison between methods (both performance and correctness),

b) investigation of a single method in different circumstances.


Develop a computational structure that can be reused with other special functions (abstract code and provide the main options), and that could also be translated to C (so: only one numerical implementation, not 100 special function numerical implementations)

Limitation: the Maple original evalf command does not accept optional arguments


The cost of numerically evaluating an Appell function



If it is a special hypergeometric case, then between 1 to 2 hypergeometric functions


Next simplest case (series/recurrence below) 3 to 4 hypergeometric functions plus adding somewhat large formulas that involve only arithmetic operations up to 20,000 times (frequently less than 100 times)


Next simplest case: the formulas themselves are power series with hypergeometric function coefficients; these cases frequently converge rapidly but may involve the numerical evaluation of up to hundreds of hypergeometric functions to get the value of a single Appell function.


Strategy for the numerical evaluation of Appell functions (or other functions ...)



The numerical evaluation flows orderly according to:

1) check whether it is a singular case

2) check whether it is a special value

3) compute the value using a series derived from a recurrence related to the underlying ODE

4) perform an sum using an infinite sum formula, checking for convergence

5) perform the numerical integration of the ODE underlying the given Appell function

6) perform a sequence of concatenated Taylor series expansions





Numerical integration of an underlying differential equation (ODEs and dsolve/numeric)


Concatenated Taylor series expansions covering the whole complex plane




Improvements in the numerical evaluation of hypergeometric functions


Evalf: an organized structure to implement the numerical evaluation of special functions in general


To be done



Download Appell_Functions.pdf

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft


I am working with numerical function with two indices, such as f[i,j], where i,j run over 1 to N, where N is n integer number.

The resluts stored in vector colummn for each value of j. i.e.

for each j=1,     G1:=Vector[column](N,i-->f(i,1))

for each j=2,     G2:=Vector[column](N,i-->f(i,2))

and so on.

The procedure is working for the first value j=1.

for the second value the following message appear

{--> enter Terminate, args = 
{--> enter Terminate, args = 
<-- exit Terminate (now in Terminate) = }
<-- exit Terminate (now at top level) = }

Hi, i try resolve this equations numerical for the boundary between the spinglass and ferromagnetic phases:






Please help me to solve the system of 1st order singular O.D.E  (see uploaded file)....New_Microsoft_Office_Word_Document.docx


I have to solve numerical trogonometric equations such as :

But, after, I would like to keep only the solution defined in a specific interval such as : [0,Pi]

1) Is there a possibility to define options with the function solve to limit the solutions belonging to a specific interval ?

2) Otherwise, may you help me to make an systematic process to choose a solution in a specific interval ?

Thank you for your help


I am trying to solve a matrix system to find the relative arrival rates of a queueing network using Gauss-Seidel.The maple commands are below:

with(Student[NumericalAnalysis]); with(LinearAlgebra);

A := Matrix([[1, -.333, -.333, -.333], [0, 1, -.333, -.333], [0, -.333, 1, -.333], [0, -.333, -.333, 1]]);

IsDefinite(A, 'query' = 'positive_semidefinite');


b := Vector([1, 1, 1, 1]);

IterativeApproximate(A, initialapprox = Vector([1, 1, 1, 1]), tolerance = 10^(-3), maxiterations = 20, stoppingcriterion = relative(infinity), method = gaussseidel);

Error, (in Student:-NumericalAnalysis:-IterativeApproximate) check that the augmented matrix has the correct dimensions

I do not understand this error as the matrix is 4x4 as shown. Can anyone see where I went wrong?


Hello Mapleprime users

I am having an issue with a numerical integration calculation. I have a large(ish) polynomial integrand and need to apply two integrations which I am using the _cuhre method for. See attached the minimum working example file.

Firstly some simplifications are done to the integrand and then a basic for loop which calculates the integration for r from 0 to 20 in steps of 0.1. The issue occurs when the loop hits r~2.5 (on my machine, Maple 2015, i5, 16GB ram). Up to that point the calculation is steady with the following stats per point calculated:

memory used=87.59MiB, alloc change=0 bytes, cpu time=3.45s, real time=3.15s, gc time=439.16ms

Then when the calculation gets to ~2.5 it just hangs and will not calculate past it. Any ideas as to what is going on here?

Any help would be appreciated with this issue.

Thank you in advance




Hi, Maple community

I am trying to integrate a 3D data numerical set. I managed to execute an array interpolation algorithm and apply it to plot a smooth surface but I do not know how to integrate it numerically. I want to be able to integrate the data not only in its domain but also in sub-domains for multiple purposes.

I've found here a topic about it and tried to apply suggestions made with no success.

Here I copy the worksheet steps:


Initial Parameters

Nx:=8;    Points in the x direction
Nz:=6;    Points in the z direction
Lx:=3;    This means x goes from 0 to 3
Lz:=2;    This means z goes from 0 to 2

DataPn:=Matrix(Nx,Nz,0):         Matrix to store values in the 3rd direction

The following step is a function that I used to simulate data that I will later on will get by measurings.

for i from 1 to Nx do
 for j from 1 to Nz do
 end do;
end do;

Setting data to use ArrayInterpolation


In this step here, I don't really understand what this argument [[a],[b]] means but ArrayInterpolation cannot work without it.

MI:=(a,b)->ArrayInterpolation([datax, dataz],Array(DataPn),[[a],[b]],'method' = 'spline')[1,1];

Plot and display the interpolation and the real function to compare the quality of the interpolation, actually it does an excellent job, at least with this data.

G1:=plot3d(MI, datax[1]..datax[Nx], dataz[1]..dataz[Nz],labels=[x,z,Pn]);

And finally I found this integration method here in mapleprimes

evalf(Int(x->evalf(Int( y->MI(x,y), 0..3,method=_d01akc)),0..2,method=_d01akc));

And the result is 1.102279199. I thought maybe I was setting wrongly the integration limits or the order in the command but in any case the result was near the analytic case which was:



In this scenario I know the result is wrong because I can compute the exact solution, but later on I won't be able to do so. Also, later I'm interested in the integration, let's say, x from 0..0.5 and later from 0.5..1 and so on (keeping z always from 0 to 2) So, the main questions are:

1. Am I setting something wrong?
2. If I delete the method option in the integral, it takes a lot of time for the program to compute any result. Why?
3. Is there any other way to integrate this? I mean, any other way to write the command to compute the integral.
4. What does the [[a],[b]] argument stands for in ArrayInterpolation?

Regards and thanks for your help.

Hello everyone,

I am trying to solve numerically int( f(t,z) , t=0..T ) = 0 , in z for a cumbersome f.

I tried z1=fsolve( int( f(t,z) , t=0..T ) = 0 , z). But then I tried int( f(t,z1) , t=0..T ) and the result is clearly not zero nor anything small.

It looks like Maple evaluates analytically the integral, and does it wrong (check this for more details) so fsolve uses the wrong equations.

Anyone knows how I can force Maple to evaluate numerically the integral at each step of the fsolve function?

Thank you!

I am solving a complicated ODE and I would like to know if there is a way for Maple to output the ODEs without doing any numerical substitutions for known parameters. Say one of my parameters, call it P, is initialized (there are many more but Ill just simplify and consider one here) to a value of 10. In order to chekc that I have coded the ODEs correctly it would help me if Maple does not substitue with the numerical value 10 for P when displayng teh ODEs, but rather keeps P, as a parameter. Is there a way to achieve this?

Hi all,

I'm trying to perform some calculations containing exponential integrals. Here is a snippet of my code:



Es := 4*10^9:


C := (q) -> q^(H+1.5):

G := (q) -> (int(q^3*C(q), q = q[0] .. q)):
P := (q) -> 1/sqrt(G(q)):

p := (xi) -> p[0]/P(xi*q[0]):
w:= (q,xi) -> 1/(Pi*(int(C(qs)*qs^3, qs = xi*q[0] .. q)))^(1/2):

U := (xi) -> (int(q*C(q)*w(q,xi)*(Int(exp(-(w(q, xi)*ps/Es)^2)/ps,ps=p(xi)..infinity)),q=xi*q[0]..q[s])):

My objective is to evaluate U for a set of discrete values of xi for further processing e.g. visualisation via plots. Neither value(U(xi)) nor evalf(U(xi)) produces a numerical result so I keep searching for solution. Does anybody have an advice how to solve U?

Regards, lassa



I wonder how is it possible to numerically evaluate two-dimensional sum, something like this:


Hello everyone !

I have a problem when I want to calculate the following multiple integration numerically:


It doesn't work. But when I replace sum(x[i],i=1..6)^2 with sum(x[i],i=1..6), it works. Is there any feasible solution to my problem ?

Thank you for reading !


My goal is to plot the integral J with respect to t and as you can see J is a piecewise function.

This is my code.


Actually it's a problem about adiabatic invariants.

If you want to know the backgroud please see this link.

I have the solution to a heat PDE, v and the error esitmate u + cos(x+t) = v


I want to plot log v(1,t) as function of log u(1,t) in maple, but I seem to get an error:

Error, (in plot) unexpected option: ln(u(1, t))


I am attaching my code below.

How to fix this problem?

Thanks in advance.

1 2 3 4 5 6 7 Page 1 of 8