Items tagged with algorithm

The 196 algorithm goes like this.  Start with an integer.  Reverse the digits.  Add the reversed number to the integer.  For most numbers, this eventually leads to a palendrome.  That is to say the number is equal to the reversed number.  I wrote a little Maple procedure to explore 196, the smallest number that will probrably never become a palendrome when put into the algorithm.


Let me know if you like my code.




In a recent conversation I explained whyLSODE was giving wrong results ( After a lot of confusions and weird infinite loops for answers, it turned out that Newton Raphson was not properly done.

Both LSODE and MEBDFI are currently incompletely implemented (only one iteration is done instead of Newton Raphson till convergence). Maplesoft should update the help files accordingly.

The post below explains how better results are obtained with method = mgear. To run the command mgear you will need Maple 6 or earlier versions. For lsode, any current version is fine.  Unfortunately Maple deprecated an algorithm that worked fine. From Maple 8, the algorithm moved to Rosenbrock methods for stiff equations. This is still not ideal.

If Maple had a working algorithm, I am hoping that Maplesoft folks would consider bringing it back in future versions. (At least with the same functionality as in Maple 6).

PLEASE NOTE, the issue is not with solving this example (Very simple). This example is chosen to show how a popular algorithm in the literature is wrongly implemented.


Here Maple's lsode is forced to take only one step and use first order back ward difference formula to integrate from 0 to 1.  LSODE mimics Eulerbackward using the options given below. The post shows that LSODE does not do Newton Raphson and just performs only iteration for nonlinear equations.



Digits := 15



eq := diff(y(t), t) = -y(t)



C := Vector[row](22, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 0, (19) = 0, (20) = 0, (21) = 0, (22) = 0})



C[9] := 1




[t = .1, y(t) = .909090909090834]



0.1e2*y1-0.1e2 = -y1





 While for linear it gave the expected result, it gives wrong results for nonlinear problems.



[t = .1, y(t) = .904837355407810]



eq := diff(y(t), t) = -y(t)^2*exp(-y(t))-10*y(t)*(1+0.1e-1*exp(y(t)))




[t = .1, y(t) = .501579294869466]



0.1e2*y1-0.1e2 = -y1^2*exp(-y1)-10*y1*(1+0.1e-1*exp(y1))






 the expected answer is correctly obtained with default tolerance as


[t = .1, y(t) = .349614721994122]


 The results obtained are worse than single iteration using jacobian.


eq2 := 0.1e2*y1-0.1e2+y1^2*exp(-y1)+10*y1*(1+0.1e-1*exp(y1))



jac := proc (y1) options operator, arrow; 20.+2*y1*exp(-y1)-y1^2*exp(-y1)+.10*exp(y1)+.10*y1*exp(y1) end proc



f := proc (y1) options operator, arrow; 0.1e2*y1-0.1e2+y1^2*exp(-y1)+10*y1*(1+0.1e-1*exp(y1)) end proc



y0 := 1



dy := -.508796088545793



ynew := .491203911454207


 Following procedures confirm that it is indeed calling the procedure only at 0 and 0.1, with backdiag giving slightly better results.

myfun:= proc(x,y) if not type(x,'numeric') or not type(evalf(y),numeric)then 'procname'(x,y);
    else lprint(`Request at x=`,x); -y^2*exp(-y(x))-10*y*(1+0.01*exp(y)); end if; end proc;

myfun := proc (x, y) if not (type(x, 'numeric') and type(evalf(y), numeric)) then ('procname')(x, y) else lprint(`Request at x=`, x); -y^2*exp(-y(x))-10*y*(1+0.1e-1*exp(y)) end if end proc




`Request at x=`, 0.

`Request at x=`, 0.

`Request at x=`, .1

`Request at x=`, .1

[x = .1, y(x) = .501579304183583]




`Request at x=`, 0.

`Request at x=`, 0.

`Request at x=`, .1

`Request at x=`, .1

[x = .1, y(x) = .497831388424072]



Download Lsodeanalysistrunc.mws


Next see how dsolve method = mgear works just fine in Maple 6 (gives the expected answer upto 3 Digits accuracy). To run this code you will need Maple 6 or earlier versions. Maple 7 has this algorithm, but I don't know to use it as it is hidden. I would like to get support from other members to get Maplesoft's attention to bring this algorithm back.

If Mdy/dt = f(y) is solved using mgear algorithm (instead of dy/dt =f ), then one can have a good DAE solver based on this (M being singular). 



myfun:= proc(x,y) if not type(x,'numeric') or not type(evalf(y),numeric)then 'procname'(x,y);
    else lprint(`Request at x=`,x); -y^2*exp(-y(x))-10*y*(1+0.01*exp(y)); end if; end proc;

myfun := proc (x, y) if not (type(x, 'numeric') and type(evalf(y), numeric)) then ('procname')(x, y) else lprint(`Request at x=`, x); -y^2*exp(-y(x))-10*y*(1+0.1e-1*exp(y)) end if end proc




`Request at x=`, 0.

`Request at x=`, .1

`Request at x=`, .1

`Request at x=`, .1

[x = .1, y(x) = .4887165263]




Download Mgearworks.mws

Greetings to all. I am writing today to share a personal story / exploration using Maple of an algorithm from the history of combinatorics. The problem here is to count the number of strings over a certain alphabet which consist of some number of letters and avoid a set of patterns (these patterns are strings as opposed to regular expressions.) This counting operation is carried out using rational generating functions that encode the number of admissible strings of length n in the coefficients of their series expansions. The modern approach to this problem uses the Goulden-Jackson method which is discussed, including a landmark Maple implementation from a paper by D. Zeilberger and J. Noonan, at the following link at (Goulden-Jackson has its own website, all the remaining software described in the following discussion is available at the MSE link.) The motivation for this work was a question at the MSE link about the number of strings over a two-letter alphabet that avoid the pattern ABBA.

As far as I know before Goulden-Jackson was invented there was the DFA-Method (Deterministic Finite Automaton also known as FSM, Finite State Machine.) My goal in this contribution was to study and implement this algorithm in order to gain insight about its features and how it influenced its powerful successor. It goes as follows for the case of a single pattern string: compute a DFA whose states represent the longest prefix of the pattern seen at the current position in the string as it is being scanned by the DFA, with the state for the complete pattern doubling as a final absorbing state, since the pattern has been seen. Translate the transitions of the DFA into a system of equations in the generating functions representing strings ending with a given maximal prefix of the pattern, very much like Markov chains. Finally solve the system of equations for the generating functions and thus obtain the sequence of values of strings of length n over the given alphabet that avoid the given pattern.

I have also implemented the DFA method for sets of patterns as opposed to just one pattern. The algorithm is the same except that the DFA does not consist of a chain with backlinks as in the case of a single pattern but a tree of prefixes with backlinks to nodes higher up in the tree. The nodes in the tree represent all prefixes that need to be tracked where obviously a common prefix between two or more patterns is shared i.e. only represented once. The DFA transitions emanating from nodes that are leaves represent absorbing states indicating that one of the patterns has been seen. We run this algorithm once it has been verified that the set of patterns does not contain pairs of patterns where one pattern is contained in another, which causes the longer pattern to be eliminated at the start. (Obviously if the shorter pattern is forbidden the so is the longer.) The number of states of the DFA here is bounded above by the sum of the lengths of the patterns with subpatterns eliminated. The uniqueness property of shared common prefixes holds for subtrees of the main tree i.e. recursively. (The DFA method also copes easily with patterns that have to occur in a certain order.)

I believe the Maple code that I provide here showcases many useful tricks and techniques and can help the reader advance in their Maple studies, which is why I am alerting you to the web link at MSE. I have deliberately aimed to keep it compatible with older versions of Maple as many of these are still in use in various places. The algorithm really showcases the power of Maple in combinatorics computing and exploits many different aspects of the software from the solution of systems of equations in rational generating functions to the implementation of data structures from computer science like trees. Did you know that Maple permits nested procedures as known to those who have met Lisp and Scheme during their studies? The program also illustrates the use of unit testing to detect newly introduced flaws in the code as it evolves in the software life cycle.

Enjoy and may your Maple skills profit from the experience!

Best regards,

Marko Riedel

The software is also available here: dfam-mult.txt

Hi everybody,

how about a little arithmetical challenge?


known quantities : Z, t, epsilon : reals

let W = X+Y*t    X, Y integers, t real (t with unlimited precision)

let W = Z + epsilon (epsilon < 1/W), Z real (Z with unlimited precision)

W and Z have same decimal expansion ( W - Integerpart[W], Z - IntegerPart[Z] ), up to a certain "rank" (determined by the precision : epsilon)

Y*t and Z have same decimal expansion, up to a certain "rank".

The problem : build an algorithm to find the integer Y, using the "usable" decimal expansion of Z.


Thanks in advance for your contributions!

Hi there, fellow primers, it's good to be back after almost 5 years! I just want to share a worksheet on Numerov's algorithm in Maple using procedures as I've recently found out that google could not find any Maple procedure that implements Numerov's algorithm to solve ODEs.   Reference.pdf 

I'm trying to implement the QR algorithm to find the Eigenvalues of the input matrix which will be forwarded to another implementation (of the SVD alg.) to find the singular values. My implementation goes as follows:

1. feeding input: A::Matrix(datatype=float) # a bidiagonal matrix
2. construct input matrix for the QR alg. of matrix A and Z (zeros of size A): C := Matrix([[Z,Transpose(A)],[A,Z]], datatype=float); # therefore C should be symmetric
3. find the eigenvalues of matrix C with an implementation of the QR alg.:

for k from 1 to 400 do
Q, R := QRDecomposition(C);
end do:

At this point, the eigenvalues of C should be placed in the diagonal of the matrix, but they're randomly placed around the diagonal, with only ~0 elements (like 2,xxx * 10^(-13)) in the diagonal.

If anyone knows how to resolve this, let the knowledge flow through. Any help will be appriciated, thanks in advance.

I have to use the optimization package. 

- The objective function is non linear,

- I have constrains and bounds,

- The constrains are not linear.


I have reading the help page on


My question are:
Can you confirm me that the only algorithm I can use is : NLPSolve with method ''sqp''?

And if I would like to use the gradient method how can I do?

I am looking for a quick implementation in Maple to approximate a solution to the travelling salesman problem in a graph of several thousand vertices randomly placed in the unit square. I've tried TravelingSaleman in the GraphTheory package; it bogs down with just a handful (say, twelve) of vertices. I also tried Bruno Buerrier's implementation of the two-opt algorithm; on my machine that took more than a day for 2048 vertices. 

Is anyone aware of a quicker implementation in Maple?

I faced a very large eigenproblem during my research. The square matrix under consideration is of size more than 2^30 times 2^30. I have tried to deal with this problem by the QR algorithm with double implicit shift (more precisely, the Francis double step QR algorithm). I'm a very beginner of programming, but I tried as follows:


A := Matrix([[7, 3, 4, -11, -9, -2], [-6, 4, -5, 7, 1, 12], [-1, -9, 2, 2, 9, 1], [-8, 0, -1, 5, 0, 8], [-4, 3, -5, 7, 2, 10], [6, 1, 4, -11, -7, -1]]):
H := HessenbergForm(A):
for p while p>2 do: 
for k from 0 to p-3 do:  
if k<3 then z:=H(k+4,k+1):   
end if: 
if abs(H(p,q))<10^(-20)*(abs(H(q,q))+abs(H(p,p))) then    H(p,q):=0: p:=p-1:q=p-1:  
elif abs(H(p-1,q-1))<10^(-20)*(abs(H(q-1,q-1))+abs(H(q,q))) then    H(p-1,q-1):=0: p:=p-2:q:=p-1:  
end if:  od:

It seemed that replacing 0 in a Hessenberg matrix by a non-zero element is not allowed. How can I remedy this?

Plus, can anyone tell me the problem of the above thing(it's not really a programming...;( ), please?

I would also appreciate it if someone let me know a better idea for a huge eigenproblem.

Thanks in advance.

In this paper we will demonstrate the importance of using simple to complex algorithms applied to complex systems in civil and mechanical engineering. In order to develop solutions that developers need to be involved in issues of advanced dynamic computer science. We show how is that with the Maple scientific program and through component-based algorithms can generate power then then be inserted into specific algorithms. Will form patterns with movements of rotation and revolution of their axes, in each case to model and analyze the curves thereof comprising. With these modelalos and curve analysis we can predict manufacturing costs, freight, inter alia estrcturas which they can be used with the correct use of Maplesoft.



(in spanish)

Lenin Araujo Castillo





Hello everbody!


uses LA=LinearAlgebra;  

local      i,k,n:= LA:-RowDimension(A),    

        x:= Vector(LA:-RowDimensions(A)),    

        p:= Vector(LA:-RowDimensions(A));  


 while  k<=n do      

       for i to n while i<>j do          

          x[i]:=1/(A[i,i])*(-add(A[i,j]*p[j],j=1..n)+b[i]);       end do;      

        if abs(x-p)<epsilon do return x; end if;    



 end do;  


end proc:

Hello! I have written a algorithm. Can you help me find errors? thank you very much. sorry, my English is not very good!

uses LA= LinearAlgebra;
local i, j, k, n:= LA:-RowDimension(A),
L:= Matrix(LA:-Dimensions(A));
for j from 2 to n do
end do;
for i from 2 to n-1 do
for j from i+1 to n do
end do;
end do;

hello eveyone! sorry, my English is not very good

I writed Neville algorithm

I want to creat a table(or a matrix) Q with



with value of X: -2,1,0,1,2

-2   1/4

-1    1/2

0    1

1    2

2     4

I want to approximate f at x=0.5 by Neville

then:   for i:=2,...,n   (that case is 5)

             for j:=2,...i




this is:

-2  1/4     0       0            0          0

-1  1/2     0.875  0           0           0

0  1         1.25   1.3475    0           0

1  2         1.5    1.4375   1.421875  0

2  4         1      1.375      1.40625  1.412109375

do you understand my mind? sorry, my English is not very good



return 2^x
for k to n+1 do
for i to n+1 do
if i<>k then
L:=L*(x-M[i])/(M[k]-M[i]); end if; end do;
end do;
I write Lagran algorithm. sorry, my english is not very good

1 2 3 4 5 6 Page 1 of 6