Main Content

Results for


syms u v
atan2alt(v,u)
ans = 
function Z = atan2alt(V,U)
% extension of atan2(V,U) into the complex plane
Z = -1i*log((U+1i*V)./sqrt(U.^2+V.^2));
% check for purely real input. if so, zero out the imaginary part.
realInputs = (imag(U) == 0) & (imag(V) == 0);
Z(realInputs) = real(Z(realInputs));
end
As I am editing this post, I see the expected symbolic display in the nice form as have grown to love. However, when I save the post, it does not display. (In fact, it shows up here in the discussions post.) This seems to be a new problem, as I have not seen that failure mode in the past.
You can see the problem in this Answer forum response of mine, where it did fail.
D.R. Kaprekar was a self taught recreational mathematician, perhaps known mostly for some numbers that bear his name.
Today, I'll focus on Kaprekar's constant (as opposed to Kaprekar numbers.)
The idea is a simple one, embodied in these 5 steps.
1. Take any 4 digit integer, reduce to its decimal digits.
2. Sort the digits in decreasing order.
3. Flip the sequence of those digits, then recompose the two sets of sorted digits into 4 digit numbers. If there were any 0 digits, they will become leading zeros on the smaller number. In this case, a leading zero is acceptable to consider a number as a 4 digit integer.
4. Subtract the two numbers, smaller from the larger. The result will always have no more than 4 decimal digits. If it is less than 1000, then presume there are leading zero digits.
5. If necessary, repeat the above operation, until the result converges to a stable result, or until you see a cycle.
Since this process is deterministic, and must always result in a new 4 digit integer, it must either terminate at either an absorbing state, or in a cycle.
For example, consider the number 6174.
7641 - 1467
ans = 6174
We get 6174 directly back. That seems rather surprising to me. But even more interesting is you will find all 4 digit numbers (excluding the pure rep-digit nmbers) will always terminate at 6174, after at most a few steps. For example, if we start with 1234
4321 - 1234
ans = 3087
8730 - 0378
ans = 8352
8532 - 2358
ans = 6174
and we see that after 3 iterations of this process, we end at 6174. Similarly, if we start with 9998, it too maps to 6174 after 5 iterations.
9998 ==> 999 ==> 8991 ==> 8082 ==> 8532 ==> 6174
Why should that happen? That is, why should 6174 always drop out in the end? Clearly, since this is a deterministic proces which always produces another 4 digit integer (Assuming we treat integers with a leading zero as 4 digit integers), we must either end in some cycle, or we must end at some absorbing state. But for all (non-pure rep-digit) starting points to end at the same place, it seems just a bit surprising.
I always like to start a problem by working on a simpler problem, and see if it gives me some intuition about the process. I'll do the same thing here, but with a pair of two digit numbers. There are 100 possible two digit numbers, since we must treat all one digit numbers as having a "tens" digit of 0.
N = (0:99)';
Next, form the Kaprekar mapping for 2 digit numbers. This is easier than you may think, since we can do it in a very few lines of code on all possible inputs.
Ndig = dec2base(N,10,2) - '0';
Nmap = sort(Ndig,2,'descend')*[10;1] - sort(Ndig,2,'ascend')*[10;1];
I'll turn it into a graph, so we can visualize what happens. It also gives me an excuse to employ a very pretty set of tools in MATLAB.
G2 = graph(N+1,Nmap+1,[],cellstr(dec2base(N,10,2)));
plot(G2)
Do you see what happens? All of the rep-digit numbers, like 11, 44, 55, etc., all map directly to 0, and they stay there, since 0 also maps into 0. We can see that in the star on the lower right.
G2cycles = cyclebasis(G2)
G2cycles = 2x1 cell array
{1x1 cell} {1x5 cell}
G2cycles{1}
ans = 1x1 cell array
{'00'}
All other numbers eventually end up in the cycle:
G2cycles{2}
ans = 1x5 cell array
{'09'} {'45'} {'27'} {'63'} {'81'}
That is
81 ==> 63 ==> 27 ==> 45 ==> 09 ==> and back to 81
looping forever.
Another way of trying to visualize what happens with 2 digit numbers is to use symbolics. Thus, if we assume any 2 digit number can be written as 10*T+U, where I'll assume T>=U, since we always sort the digits first
syms T U
(10*T + U) - (10*U+T)
ans = 
So after one iteration for 2 digit numbers, the result maps ALWAYS to a new 2 digit number that is divisible by 9. And there are only 10 such 2 digit numbers that are divisible by 9. So the 2-digit case must resolve itself rather quickly.
What happens when we move to 3 digit numbers? Note that for any 3 digit number abc (without loss of generality, assume a >= b >= c) it almost looks like it reduces to the 2 digit probem, aince we have abc - cba. The middle digit will always cancel itself in the subtraction operation. Does that mean we should expect a cycle at the end, as happens with 2 digit numbers? A simple modification to our previous code will tell us the answer.
N = (0:999)';
Ndig = dec2base(N,10,3) - '0';
Nmap = sort(Ndig,2,'descend')*[100;10;1] - sort(Ndig,2,'ascend')*[100;10;1];
G3 = graph(N+1,Nmap+1,[],cellstr(dec2base(N,10,2)));
plot(G3)
This one is more difficult to visualize, since there are 1000 nodes in the graph. However, we can clearly see two disjoint groups.
We can use cyclebasis to tell us the complete story again.
G3cycles = cyclebasis(G3)
G3cycles = 2x1 cell array
{1x1 cell} {1x1 cell}
G3cycles{:}
ans = 1x1 cell array
{'000'}
ans = 1x1 cell array
{'495'}
And we see that all 3 digit numbers must either terminate at 000, or 495. For example, if we start with 181, we would see:
811 - 118
ans = 693
963 - 369
ans = 594
954 - 459
ans = 495
It will terminate there, forever trapped at 495. And cyclebasis tells us there are no other cycles besides the boring one at 000.
What is the maximum length of any such path to get to 495?
D3 = distances(G3,496) % Remember, MATLAB uses an index origin of 1
D3 = 1x1000
Inf 6 5 4 3 1 2 3 4 5 6 6 5 4 3 1 2 3 4 5 5 5 5 4 3 1 2 3 4 5
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
D3(isinf(D3)) = -inf; % some nodes can never reach 495, so they have an infinite distance
plot(D3)
The maximum number of steps to get to 495 is 6 steps.
find(D3 == 6) - 1
ans = 1x54
1 10 11 100 101 110 112 121 122 211 212 221 223 232 233 322 323 332 334 343 344 433 434 443 445 454 455 544 545 554
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
So the 3 digit number 100 required 6 iterations to eventually reach 495.
shortestpath(G3,101,496) - 1
ans = 1x7
100 99 891 792 693 594 495
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
I think I've rather exhausted the 3 digit case. It is time now to move to the 4 digit problem, but we've already done all the hard work. The same scheme will apply to compute a graph. And the graph theory tools do all the hard work for us.
N = (0:9999)';
Ndig = dec2base(N,10,4) - '0';
Nmap = sort(Ndig,2,'descend')*[1000;100;10;1] - sort(Ndig,2,'ascend')*[1000;100;10;1];
G4 = graph(N+1,Nmap+1,[],cellstr(dec2base(N,10,2)));
plot(G4)
cyclebasis(G4)
ans = 2x1 cell array
{1x1 cell} {1x1 cell}
ans{:}
ans = 1x1 cell array
{'0000'}
ans = 1x1 cell array
{'6174'}
And here we see the behavior, with one stable final point, 6174 as the only non-zero ending state. There are no circular cycles as we had for the 2-digit case.
How many iterations were necessary at most before termination?
D4 = distances(G4,6175);
D4(isinf(D4)) = -inf;
plot(D4)
The plot tells the story here. The maximum number of iterations before termination is 7 for the 4 digit case.
find(D4 == 7,1,'last') - 1
ans = 9985
shortestpath(G4,9986,6175) - 1
ans = 1x8
9985 4086 8172 7443 3996 6264 4176 6174
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
Can you go further? Are there 5 or 6 digit Kaprekar constants? Sadly, I have read that for more than 4 digits, things break down a bit, there is no 5 digit (or higher) Kaprekar constant.
We can verify that fact, at least for 5 digit numbers.
N = (0:99999)';
Ndig = dec2base(N,10,5) - '0';
Nmap = sort(Ndig,2,'descend')*[10000;1000;100;10;1] - sort(Ndig,2,'ascend')*[10000;1000;100;10;1];
G5 = graph(N+1,Nmap+1,[],cellstr(dec2base(N,10,2)));
plot(G5)
cyclebasis(G5)
ans = 4x1 cell array
{1x1 cell} {1x2 cell} {1x4 cell} {1x4 cell}
ans{:}
ans = 1x1 cell array
{'00000'}
ans = 1x2 cell array
{'53955'} {'59994'}
ans = 1x4 cell array
{'61974'} {'63954'} {'75933'} {'82962'}
ans = 1x4 cell array
{'62964'} {'71973'} {'83952'} {'74943'}
The result here are 4 disjoint cycles. Of course the rep-digit cycle must always be on its own, but the other three cycles are also fully disjoint, and are of respective length 2, 4, and 4.
This stems purely from some play on my part. Suppose I asked you to work with the sequence formed as 2*n*F_n + 1, where F_n is the n'th Fibonacci number? Part of me would not be surprised to find there is nothing simple we could do. But, then it costs nothing to try, to see where MATLAB can take me in an explorative sense.
n = sym(0:100).';
Fn = fibonacci(n);
Sn = 2*n.*Fn + 1;
Sn(1:10) % A few elements
ans = 
For kicks, I tried asking ChatGPT. Giving it nothing more than the first 20 members of thse sequence as integers, it decided this is a Perrin sequence, and gave me a recurrence relation, but one that is in fact incorrect. Good effort from the Ai, but a fail in the end.
Is there anything I can do? Try null! (Look carefully at the array generated by Toeplitz. It is at least a pretty way to generate the matrix I needed.)
X = toeplitz(Sn,[1,zeros(1,4)]);
rank(X(5:end,:))
ans = 5
Hmm. So there is no linear combination of those columns that yields all zeros, since the resulting matrix was full rank.
X = toeplitz(Sn,[1,zeros(1,5)]);
rank(X(6:end,:))
ans = 5
But if I take it one step further, we see the above matrix is now rank deficient. What does that tell me? It says there is some simple linear combination of the columns of X(6:end,:) that always yields zero. The previous test tells me there is no shorter constant coefficient recurrence releation, using fewer terms.
null(X(6:end,:))
ans = 
Let me explain what those coefficients tell me. In fact, they yield a very nice recurrence relation for the sequence S_n, not unlike the original Fibonacci sequence it was based upon.
S(n+1) = 3*S(n) - S_(n-1) - 3*S(n-2) + S(n-3) + S(n-4)
where the first 5 members of that sequence are given as [1 3 5 13 25]. So a 6 term linear constant coefficient recurrence relation. If it reminds you of the generating relation for the Fibonacci sequence, that is good, because it should. (Remember I started the sequence at n==0, IF you decide to test it out.) We can test it out, like this:
SfunM = memoize(@(N) Sfun(N));
SfunM(25)
ans = 3751251
2*25*fibonacci(sym(25)) + 1
ans = 
3751251
And indeed, it works as expected.
function Sn = Sfun(n)
switch n
case 0
Sn = 1;
case 1
Sn = 3;
case 2
Sn = 5;
case 3
Sn = 13;
case 4
Sn = 25;
otherwise
Sn = Sfun(n-5) + Sfun(n-4) - 3*Sfun(n-3) - Sfun(n-2) +3*Sfun(n-1);
end
end
A beauty of this, is I started from nothing but a sequence of integers, derived from an expression where I had no rational expectation of finding a formula, and out drops something pretty. I might call this explorational mathematics.
The next step of course is to go in the other direction. That is, given the derived recurrence relation, if I substitute the formula for S_n in terms of the Fibonacci numbers, can I prove it is valid in general? (Yes.) After all, without some proof, it may fail for n larger than 100. (I'm not sure how much I can cram into a single discussion, so I'll stop at this point for now. If I see interest in the ideas here, I can proceed further. For example, what was I doing with that sequence in the first place? And of course, can I prove the relation is valid? Can I do so using MATLAB?)
(I'll be honest, starting from scratch, I'm not sure it would have been obvious to find that relation, so null was hugely useful here.)
I saw this post on Answers.
I was impressed at the capability of the AI, as I have been at other times when I posed a question to it, at least some of the time. So much so that I wondered...
What if the AI were automatically applied to EVERY question on Answers? Would that be a good or bad thing? For example, suppose the AI automatically offers an answer to every question as soon as it gets posted? Of course, users would still be allowed to post their own, possibly better answers. But would it tend to disincentivise individuals from ansering questions?
Perhaps as bad, would it push Answers into the mode of a homework solving forum? Since if every homework question gets a possibly pretty good automatic AI generated solution, then every student will just post all HW questions, and the forum would quickly become overwhelmed.
I suppose one idea could be to set up the AI to post an answer to all un-answered questions that are at least one month old. Then students would not gain by posting their homework.
I've now seen linear programming questions pop up on Answers recently, with some common failure modes for linprog that people seem not to understand.
One basic failure mode is an infeasible problem. What does this mean, and can it be resolved?
The most common failure mode seems to be a unbounded problem. What does this mean? How can it be avoided/solved/fixed? Is there some direction I can move where the objective obviously grows without bounds towards +/- inf?
Finally, I also see questions where someone wants the tool to produce all possible solutions.
A truly good exposition about linear programming would probably result in a complete course on the subject, and Aswers is limited in how much I can write (plus I'll only have a finite amount of energy to keep writing.) I'll try to answer each sub-question as separate answers, but if someone else would like to offer their own take, feel free to do so as an answer, since it has been many years for me since I learned linear programming.

I saw this problem online recently.

Passenger distribution on a train

Not a terribly difficult problem to solve. But it was mildly interesting to find a solution using MATLAB. Perhaps just as interesting is the post analysis of the problem to understand what is happening, and why any unique solution exists at all for one specific car.

The question is, we have a passenger train with 11 cars in it. Feel free to number them 1 through 11. We know that 381 passengers boarded the train. Every passenger is in one of the cars, but all we know is there are exactly 99 passengers in every set of three consecutive cars. Now the question becomes, how many passengers are in car number 9?

One might say at first this is impossible to know. Surely there are many ways the passengers may be arranged, but is that true? Could it be impossible to solve?

First, before we go any further, a few tests seem important. Logically, we might think to distribute 33 passengers in every car. Would that work? So we would have a passenger distribution of

X1 = repmat(33,11,1)
X1 =
    33
    33
    33
    33
    33
    33
    33
    33
    33
    33
    33

While that satisfies the requirement of every 3 consecutive cars having 99 passengers, it fails the total count requirement, since we can see the sum of all passengers would be 363. This yields too few total passengers, with only a combined load of 363 passengers, and we need 381.

At the other end of the spectrum is another extreme case. We might have this distribution:

X2 = zeros(11,1);
X2(1:3:11) = 99
X2 =
    99
     0
     0
    99
     0
     0
    99
     0
     0
    99
     0

Again, it meets the requirement that the sum of passengers in any 3 consecutiuve cars will be 99. But that case yields too many total passengers at 396. Somewhere in the middle must/might/may be a solution, right? At least it is good to see that we can have more or less than 381 total passengers. But how can we find a solution using MATLAB?

There is one other problem with the X2 attempt at a solution, in that had I chosen a different first car to place the 99 passengers, we need not have a unique result in car number 9.

X3 = zeros(11,1);
X3(2:3:11) = 99;
X4 = zeros(11,1);
X4(3:3:11) = 99;

Each of those schemes would put 99 passengers in every set of 3 consecutive cars.

[X2,X3,X4]
ans =
    99     0     0
     0    99     0
     0     0    99
    99     0     0
     0    99     0
     0     0    99
    99     0     0
     0    99     0
     0     0    99
    99     0     0
     0    99     0

But car number 9 would have very different numbers of passengers, depending on the choice made, either 0 or 99 passengers.

An obvious solution is to look for a code that can solve such a problem for us. INTLINPROG stands out as the perfect tool, as this is a linear problem, with everything being in the form of a sum. The unknowns will be how many passengers are sitting in each car. There are 11 cars. So there are 11 unknowns. The bounds are simple.

lb = zeros(1,11); % There cannot be less than zero passengers in any car.
ub = repmat(99,1,11); % since the sum of any three consecutive cars is 99, we cannot have more than 99 people in any one car.

All of the unknown car counts must be integer. That is, we cannot have a fractional number of people in a car, unless this is part of an Agatha Christie murder mystery.

intcon = 1:11;

What constraints apply? First, the sum of all passengers on the train must be 381.

As well, we know that in every 3 consecutive cars, the sum must be 99. Both constraints will take the form of exact linear equality constraints. We can encode all of that into the matrix Aeq, and the vector Beq.

Aeq = [ones(1,11);triu(tril(ones(9,11),2))]; % A tricky way to create the matrix Aeq
Beq = [381;repmat(99,9,1)];

If X is a potential solution that satisfies the bound constraints, it must satisfy the matrix equation Aeq*X==Beq.

We can see for example, the potential solutions I posed above as X1,...X4, all fail to satisfy the requirement on the total number of passengers, since while the sums for consecutive cars are correct, the total sum is not.

Aeq*[X1,X2,X3,X4]
ans =
   363   396   396   297
    99    99    99    99
    99    99    99    99
    99    99    99    99
    99    99    99    99
    99    99    99    99
    99    99    99    99
    99    99    99    99
    99    99    99    99
    99    99    99    99

Finally, there are no linear inequality constraints.

A = [];
B = [];

At this point, you might be wondering how we can formulate this as a linear programming problem at all. What would be the objective function? What could we hope to minimize? As it turns out, linear programming tools are pretty simple in that respect. We could pose just about any objective we want. For example, this is sufficient:

f = ones(1,11);

Effectively, we are just using intlinprog to see if a FEASIBLE integer solution exists that satisfies all of the bounds, as well as the equality constraints. This is why the objective can be the same as one of the equality constraints. Once INTLINPROG finds any solution, it will be done.

And now we can throw the problem into INTLINPROG, hoping something intelligent falls out.

[X,~,EXITFLAG] = intlinprog(f,intcon,A,B,Aeq,Beq,lb,ub)
P:                Optimal objective value is 381.000000.                                           
Optimal solution found.
Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value,
options.AbsoluteGapTolerance = 0 (the default value). The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05 (the default value).
X =
     0
    84
    15
     0
    84
    15
     0
    84
    15
     0
    84
EXITFLAG =
     1

intlinprog has found a solution,

isequal(Aeq*X,Beq)
ans =
  logical
   1

The solution satisfies all of the constraints. Effectively, we see the repeating sub-sequence [0 84 15] in consecutive cars. And of course, as long as we repeat that sequence, it does satisfy all requirements. How many people are sitting in car number 9? 15 people.

Symmetry would suggest that car number 3 must have the same number of people, since we could as easily have numbered the cars starting from either end. And of course, X(3) was also 15.

Thankfully our intuition works there. It would seem we are done now. Or are we?

For some of you, you might be wondering if any other solutions can possibly exist. And some of you might be wondering if any of those solutions can have some other number of passengers than exactly 15 in cars number 3 and 9.

NULL is the MATLAB function to come to the rescue here.

This is essentially a linear algebra question. We wish to know the solutions of the problem Aeq*x==Beq. Here, Aeq is a 10x11 matrix, so it has rank at most 10. That means there is a vector Y, such that Aeq*Y == 0.

Y = double(null(sym(Aeq)))
Y =
    -1
     1
     0
    -1
     1
     0
    -1
     1
     0
    -1
     1

What does that tell us? If we have some particular solution (X) to the non-homogeneous problem Aeq*X==Beq, then the set of all possible solutions will be of the general form

syms t
X + t*Y
ans =
    -t
t + 84
    15
    -t
t + 84
    15
    -t
t + 84
    15
    -t
t + 84

You may see that this generates all solutions to the general problem. We can see a few of them in this array:

[X-1*Y, X - 2*Y, X - 3*Y, X - 84*Y]
ans =
     1     2     3    84
    83    82    81     0
    15    15    15    15
     1     2     3    84
    83    82    81     0
    15    15    15    15
     1     2     3    84
    83    82    81     0
    15    15    15    15
     1     2     3    84
    83    82    81     0

We see now there are 85 such possible integer solutions, all of the form X-k*Y, where k can be any positive integer from 0 to 84 inclusive. INTLINPROG found one of them. But as importantly, do you see that since the elements

Y([3 6 9])
ans =
     0
     0
     0

are all identically zero, that those elements in the solution can never change? Those cars must always contain exactly 15 passengers for all of the constraints to be satisfied. I'll be honest, it is not at all obvious as to why it works out that way, at least not initially in my eyes. That leaves my intuition wanting, just a bit.

How might we analyze this problem in a different way? Perhaps a different approach would yield a more satisfying solution. Suppose we chose a passenger partitioning that is strictly repetitive? For example, choose three non-negative integers u,v,w, such that u+v+w=99.

Now, fill the cars using the sequence

syms u v w
X = [u v w u v w u v w u v];

Surely you would agree that any subset of 3 consecutive cars adds to 99, as long as u+v+w=99. But then the sum of all 11 cars in that sequence must be 4*u+4v+3*w. And this leaves us with now two equations in three unknowns. We have

EQ1 = sum(X) == 381;
EQ2 = u + v + w == 99;

The rest is easy now, as we can do

EQ1 - 3*EQ2
ans =
u + v == 84

So if a solution in this form exists, we can see that u+v=84, and therefore w=99-u-v=15. (Remember that w was the number of passengers in car number 9, but also in cars numbered 3 and 6.) Any combination of non-negative integers that sums to 84 will work for u and v though.

This constructive approach does not insure it is the ONLY solution, since I built it from the sequence in the vector X. Perhaps a solution exists that is not simply repetitive as I created it. In fact, the previous analysis using null told us the whole story.