## Improved Josephus Algorithm in Maple. Need help....

Hello, everyone. I have a group project where we have to explain the Josephus problem and use Maple to solve the problem. I am trying to solve the problem in multiple ways (because why not), but I am struggling with my third procedure. I understand the logic behind it and how its supposed to achieve O(k*logn), but the code that I wrote for it doesn't seem to produce the correct result.

```JosephusImproved := proc (n, k)
local count, result:
if n = 1 then
return 0:
elif 1 < n < k then
return JosephusImproved(n - 1, k) + k + 1 mod n:
else
count := floor(n / k):
result := JosephusImproved(n - count, k):
result := result - n mod k:
if result < 0 then
result := result + n
else
result := floor(result /(k - 1)):
return result:
end if:
end if:
end proc:```

Note: The regular recursive expression [Josephus(n - 1, k) + k + 1 mod n] has a "+ 1" since that was the only way I could make Maple do the calculation correctly. Proven with a Cyclic procedure I already made.
Note 2: I am using Maple 2016 and 2D Math.

I would like some insight as to how I could fix this so that it works, just like the regular recursive procedure and cyclic list that I have.

Cheers.

## Random walk problem...

I have made an algorithm for producing random walks (only possible to walk one step to either direction except downwards(EDIT: It is possible to go downwards, but it has to be when making a turn to the right or left). The walks are determined by the rand function:

R3:=rand(1..3): (1: go straight on; 2: turn right; 3: turn left)
M:=15; N:=1500;

Randwalk3:=proc(R3)
local i,j,r,X,Y,L;
for j from 1 to M do
X[0,j]:=0;                                # Initialization
Y[0,j]:=0;
X[1,j]:=1;                                # The first step should still be taken to the point (1,0)
Y[1,j]:=0;
for i from 2 to N do
r:=R3();
if r=1 then X[i,j]:=2*X[i-1,j]-X[i-2,j]; Y[i,j]:=2*Y[i-1,j]-Y[i-2,j];                   # go straight on
elif r=2 then X[i,j]:=X[i-1,j]+Y[i-1,j]-Y[i-2,j]; Y[i,j]:=Y[i-1,j]-X[i-1,j]+X[i-2,j];    # turn right
else X[i,j]:=X[i-1,j]-Y[i-1,j]+Y[i-2,j]; Y[i,j]:=Y[i-1,j]+X[i-1,j]-X[i-2,j];             # turn left
end if;
if (X[i,j]=X[j,j] and Y[i,j]=Y[j,j]) then L[j]:=i; break; end if; (This is wrong)
end do;
end do:
return [X,Y,L];
end proc:

The question from is like this:
Modify the algorithm such that it stops at r[i] if r[i] = r[j] for any 0 <= j <= i-2. r=(xi,yi). M is the number of random walks and N is the number of steps. The length of the paths should be stored in L[m] (m=1..M). How do I implement the if-test correctly?

## Please bring back dsolve numeric method = mgear

Maple

In a recent conversation I explained whyLSODE was giving wrong results (http://www.mapleprimes.com/questions/210948-Can-We-Trust-Maple#comment230167). 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.

 > restart;
 > Digits:=15;
 (1)
 > eq:=diff(y(t),t)=-y(t);
 (2)
 > C:=array([0\$22]);
 (3)
 > C[9]:=1;
 (4)
 > sol:=dsolve({eq,y(0)=1},type=numeric,method=lsode[backfull],ctrl=C,initstep=0.1,minstep=0.1,abserr=1,relerr=1):
 > sol(0.1);
 (5)
 > subs(diff(y(t),t)=(y1-1)/0.1,y(t)=y1,eq);
 (6)
 > fsolve(%,y1=0.5);
 (7)

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

 > sol1:=dsolve({eq,y(0)=1},type=numeric):
 > sol1(0.1);
 (8)
 > eq:=diff(y(t),t)=-y(t)^2*exp(-y(t))-10*y(t)*(1+0.01*exp(y(t)));
 (9)
 > sol:=dsolve({eq,y(0)=1},type=numeric,method=lsode[backfull],ctrl=C,initstep=0.1,minstep=0.1,abserr=1,relerr=1):
 > sol(0.1);
 (10)
 > subs(diff(y(t),t)=(y1-1)/0.1,y(t)=y1,eq);
 (11)
 > fsolve(%,y1=1);
 (12)
 > sol1:=dsolve({eq,y(0)=1},type=numeric):
the expected answer is correctly obtained with default tolerance as
 > sol1(0.1);
 (13)

The results obtained are worse than single iteration using jacobian.

 > eq2:=(lhs-rhs)(subs(diff(y(t),t)=(y1-1)/0.1,y(t)=y1,eq));
 (14)
 > jac:=unapply(diff(eq2,y1),y1);
 (15)
 > f:=unapply(eq2,y1);
 (16)
 > y0:=1;
 (17)
 > dy:=-evalf(f(y0)/jac(y0));
 (18)
 > ynew:=y0+dy;
 (19)

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;
 (20)
 > sol1:=dsolve({diff(y(x),x)=myfun(x,y(x)),y(0)=1},numeric,method=lsode[backfull],ctrl=C,initstep=0.1,minstep=0.1,abserr=1,relerr=1,known={myfun}):
 > sol1(0.1);
 `Request at x=`, 0. `Request at x=`, 0. `Request at x=`, .1 `Request at x=`, .1 (21)
 > sol2:=dsolve({diff(y(x),x)=myfun(x,y(x)),y(0)=1},numeric,method=lsode[backdiag],ctrl=C,initstep=0.1,minstep=0.1,abserr=1,relerr=1,known={myfun}):
 > sol2(0.1);
 `Request at x=`, 0. `Request at x=`, 0. `Request at x=`, .1 `Request at x=`, .1 (22)

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).

 > restart;
 > 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;
 (1)
 > sol2:=dsolve({diff(y(x),x)=myfun(x,y(x)),y(0)=1},{y(x)},numeric,method=mgear[mstepnum],stepsize=0.1,minstep=0.1,errorper=1):
 > sol2(0.1);
 `Request at x=`, 0. `Request at x=`, .1 `Request at x=`, .1 `Request at x=`, .1 (2)
 >

## From combinatorics: an algorithm to enumerate...

by: Maple

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 math.stackexchange.com (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

## A little arithmetical challenge...

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.

## Numerov algorithm in Maple

by: Maple 18

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.   numerov.mw   Reference.pdf

## Problem with QR algorithm...

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);
C:=R.Q;
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.

## Approximate solution of Traveling Salesman Problem...

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?

## Eigenvalues of a huge matrix...

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):
p:=6:
for p while p>2 do:
q:=p-1:
s:=H(q,q)+H(p,p):
t:=H(q,q)*H(p,p)-H(q,p)*H(p,q):
x:=(H(1,1))^(2)+H(1,2)*H(2,1)-s*H(1,1)+t:
y:=H(2,1)*(H(1,1)+H(2,2)-s):
z:=H(2,1)*H(3,2):
for k from 0 to p-3 do:
V:=Vector([x,y,z]):
P:=Transpose(HouseholderMatrix(1/(Norm(V+exp(argument(V(1))*I)*Norm(V,2)*Vector(3,shape=unit[1]),2))*(V+exp(argument(V(1))*I)*Norm(V,2)*Vector(3,shape=unit[1])))):
r:=max(1,k):
H[k+1..k+3,r..6]:=MatrixMatrixMultiply(Transpose(P),SubMatrix(H,[k+1..k+3],[r..6])):
r:=min(k+4,6):
H[1..r,k+1..k+3]:=MatrixMatrixMultiply(SubMatrix(H,[1..r],[k+1..k+3]),P):
x:=H(k+2,k+1):
y:=H(k+3,k+1):
if k<3 then z:=H(k+4,k+1):
end if:
od:
P:=GivensRotationMatrix(Vector([x,y]),1,2):
H[q..p,p-2..6]:=MatrixMatrixMultiply(Transpose(P),SubMatrix(H,[q..p],[p-2..6])):
H[1..p,p-1,p]:=MatrixMatrixMultiply(SubMatrix(H,[1..p],[p-1,p]),P):
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.

## What is wrong with algorithm?...

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

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

## Can you help me with the program...

hello eveyone! sorry, my English is not very good

I writed Neville algorithm

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

example:

f:=X->2^X;

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

Q[i,j]=(x-X[i-j])*Q[i,j-1]-(x-X[i])*Q[i-1,j-1])/(X[i]-X[i-j])

output(Q)

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

Regards

sunflower

## What is wrong with my code?...

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

## Division of estate...

Wonder if this can be accomplished in Maple.

so I have a list of 100 items labeled {1..100} of various value {\$100, \$160, \$220, ......  , }

the task is to distribute these items among 3 people A,B,C so they get an approximately equal share.

Adding the values and dividing by 3 gives the dollar total to aim for.

This post has C.Love procedure for evenly sized groups

http://www.mapleprimes.com/questions/200480-Product-Grouping

but what i want is a method for different sized groups. ie 25 items for A, 35 for B and 40 for C (user defined).

additionally there is a fixed constraint: A has been bequeathed items 1,4,8; B items 2 and 20; C item 50.

restart:
S:= {3, 4, 5, 6, 8, 9, 28, 30, 35}:
SL:= [A,B,C,D,E,F,G,H,I]:
assign(Labels ~ (S) =~ SL); #Create remember table.
AllP:= [seq(P, P= Iterator:-SetPartitions(S, [[3,3]], compile= false))]:
lnp:= evalf(ln((`*`(S[]))^(1/3))):

Var:= proc(P::({list,set}(set)))
local r:= evalf(`+`(map(b-> abs(ln(`*`(b[]))-lnp), P)[]));
end proc:

Min:= proc(S::{list,set}, P::procedure)
local M:= infinity, X:= (), x, v;
for x in S do
v:= P(x);
if v < M then  M:= v;  X:= x  end if
end do;
X
end proc:

ans:= Min(AllP, Var);
[{3, 9, 35}, {4, 8, 28}, {5, 6, 30}]
subsindets(ans, posint, Labels);
[{I, A, F}, {B, E, G}, {C, D, H}]

## Generation of math apps in dynamic systems engineering...

Maple 2015

ABSTRACT. In this paper we demonstrate how the simulation of dynamic systems engineering has been implemented with graphics software algorithms using maple and MapleSim. Today, many of our researchers the computational modeling performed by inserting a piece of code from static work; with these packages we have implemented through the automation components of kinematics and dynamics of solids simple to complex.

It is very important to note that once developed equations study; recently we can move to the simulation; to thereby start the physical construction of the system. We will use mathematical and computational methods using the embedded buttons which lie in the dynamics leaves and viewing platform cloud of Maplesoft and power MapleNet for online evaluation of specialists in the area. Finally they will see some work done; which integrate various mechanical and computational concepts implemented for companies in real time and pattern of credibility.

Selasi_2015.pdf

(in spanish)

Lenin Araujo Castillo

 1 2 3 4 5 Page 1 of 5
﻿