Generalised solution to a n-order system of equations

I want to solve (for the first term) of a system of n-1 equations, in terms of n. How should I do it?

Essentially I want to solve for a_1 in terms of n:

(This holds for all k=1,...,n-1).

For example, with n=4, I have these 3 equations I would like to solve:

.

Once can check that satifies this set of equations. In my example, I only want .

As a more general example, this is the set of (n-1) equations:

I can write the code to solve for specific n, (like n=4 above). However, I would like to have a general formula in terms of n.

How do I code this in MatLab? I always meet a lot of errors complaining about the general formula. Thank you.

[Torsten]

Please show your code so far.

[John D'Errico]

While you said "essentially", you are not giving any essential information. Even as a mathematician myself, I do not follow what you are asking.

For example, it seems like you are talking about a system of n equations (or maybe only n-1 equations), with unknowns a_1,..., a_n. The sum goes from 1 to the smaller of n-1 or 2*k, whichever is smaller. But k seems to vary from 1 to n. Confusing even so far.

Please give a simple example. Perhaps with n = 4. Write it out CLEARLY, so someone can help you. Don't write it in terms of sums. WRITE OUT THE ACTUAL EQUATIONS YOU NEED TO SOLVE! Give someone an example, a pattern we can follow and understand.

If you cannot do that, explaining what you need to someone else, how can you or anyone possibly write code to do what has no explanation? (I would suggest this is your problem. Perhaps you do understand what you need to do, we do not know that is true or not. But I think you have only an expression there, that you are then hoping to write directly into MATLAB code using sums, etc., and have a magical solution appear.)

Do this as a comment on your question, or by editing your question to add that information.

[David Goodmanson]

Hi Aden, this seems to be defined adequately (I'm assuming that 2k on the right hand side is part of the summation), is it a homework problem?

[Aden]

Hi, Sorry for being unclear; 2k is not part of the summation. This is not a homework problem, it's part of my research. Thank you!

[Torsten]

@Aden

You could produce a graph a(1) over n, but I don't see how to get a(1) as a symbolic function of n.

[Aden]

@Torsten Here is my code, that provies a range of answers from a to b, inclusive.

% Example usage:
rangeSolutions(50,60);

n = 50:
  a1 = 24658466561372956606072219961585730241/200607444499636329906942199843718400

n = 51:
  a1 = 12095999817020908880806839172705842001/96381539428995418649814967947304320

n = 52:
  a1 = 296876662512585342903233000603819451373/2320227363430068965267573684590128000

n = 53:
  a1 = 145731168343163320419788860702885530923/1116316881358922544654288513268944000

n = 54:
  a1 = 131520038262372525763519193012546565211519/988946436362361212004404709873169728000

n = 55:
  a1 = 581429164068787472552498392939187591167489/4288711995721745855579202703142204736000

n = 56:
  a1 = 6172537074853875344545662543105451222283579/44719649242486541205586769842909654579200

n = 57:
  a1 = 30337972180019584268050518807218304758008109/215748762540829378179114736330044825456000

n = 58:
  a1 = 50635315839505187203721653607779887624851161/353936768935701591008964408883034189952000

n = 59:
  a1 = 821748863224299748123697346539124438846061261/5642413690350518644235221494358721275968000

n = 60:
  a1 = 2767201210279121773176706567854579026542586813041/18685182779533599510203630089217208636299136000

function rangeSolutions(a, b)
    % Loop from A to B
    for n = a:b
        aSym = sym('a', [n-1, 1]);  
        eqns = sym(zeros(n-1, 1));
        
        % Create the equations
        for k = 1:n-1
            eqns(k) = 2*k*aSym(k) == sum(aSym(1:min(n-1, 2*k))) + 2*k;
        end
        
        % Solve
        sol = solve(eqns, aSym);
        
        % Check if the solution is empty
        if isempty(sol)
            a1Solution = 'No solution';
        else
            % Extract a1
            a1Solution = char(sol.a1); % Convert to string
        end
        
        % Print the solution "n = (n): a1 = (a1)"
        fprintf('n = %d:\n', n);
        fprintf('  a1 = %s\n', a1Solution);
        fprintf('\n');
    end
end

[John D'Errico]

Thank you for your clarification. It was not at all clear what you meant at first, but I see your question now. It is actually an interesting one, IMHO, that gets into the linear algebra of how you solve such a question.

I'll post an answer now, first, in terms of how I would solve it for specific n, and then look to see how we might generalize it for any n.

Answer 1

[David Goodmanson]

HI Aden,

Here is a method that produces, within 4.5 %, the experimental value a_1 = 2.518. The equation is (for simplicty n-1 is replaced by n since it is the same problem)

a_k = 1 + (1/2k) Sum{j = 1,min(n,2k)} a_j     k = 1...n    (1)

The value of a_k depends on all of the a's up to a_2k. Suppose that there are no restrictions set by n, so that the upper sum limit of 2k is always allowed. Then for any constant C the solution

a_k = nC -2(k-1)                                           (2)

is exact. So for smaller values of k, that would appear to be the solution. That can't be the whole story, though, because the boundary conditions up at k = n seem to propagate down and have an effect. The lower half of the plot ends up being almost, but not quite linear.

For k >n/2 the sum index is fixed from 1 to n and the sum is a constant. Then from (1),

a_k = 1 + D/k                         exact                (3)  

for some constant D is exact for all k > n/2.

The plot of (a_k) /n for n=1000 is the blue line here. Experimentally it is almost linear all the way up to k = n/2 and then has the 1/k behavior of (3). Its slope in the almost-linear region is close to [ -2 / 1000 ] as expected.

The value of D can be determined experimentally from the upper half of the blue line, where D = k*(a_k -1) should be a constant for n/2<k<n, which it is.

The approximate calculation done below gives the red and yellow lines in the plot. For this calculation, from (2)

a_1     = nC                                      (4)
a_(n/2) = nC - 2(n/2 -1) ~~ n(C-1) 

where we start retaining only the highest powers of n. Also from (3)

a_(n/2) = 1 + D/(n/2) ~~ D/(n/2) 
a_n    =  1 + D/n     ~~ D/n

Dropping the 1 is justified since it's already known that a(n/2) ~~ n(C-1). Equating the two values of a_(n/2) gives

D = (n^2/2)(C-1)

Now C can be determined by comparing both sides of (1) for k = n. From (3) the left hand side is

 a_n = ~~ D/n = (n/2)(C-1)                          (5)

On the right hand side, the sum of the trapezoidal portion (using (4) and not including the 1/2n factor) is

 ~~ (1/2)(nC +n(C-1)) (n/2) = (n^2/4)(2C-1) 

and the sum for the 1/k portion is by approximate integration

 D Sum{k=n/2,n} (1/k) ~~ D log(2) ~~ (n^2/2) (C-1) log(2) 

Equating both sides and now including the 1/2n factor

 (n/2)(C-1) = (1/2n) ( (n^2/4)(2C-1) + (n^2/2) (C-1) log(2)

solve for C:

 C = (3/2-log(2))/(1-log(2)) = 2.6294

and from (4),

a_1 /n = C

which is too large by about 4.5%.

[Aden]

I see. To clarify, is it correct that C and D are constants over k, but are not constants over n?

I have a few questions. For simplicity now, lets assume n is even.

1. Is what you're implying, that for every n, we have ?
2. If above is true, I tested that in Wolfram Alpha for the case , (and I set k=x and n=1000 in the picture), I get the picture attached, in which the terms of x do not cancel. (I also tested it by hand, it doesn't cancel too.) Does this imply the formula is wrong?
3. In conclusion, I believe the formula holds for and possibly some small k, but it does not hold for many values of .

Edit: I don't think that the graph is linear even for small k values. If possible, could you provide a proof?

I ran code for n=28 and I obtain the set of values:

[70.17290511 68.17290511 66.16698578 64.17882444 62.15348399 60.15681025

58.24438639 56.18429445 54.07618717 52.02805718 50.08482326 48.26205992

46.56145514 44.97823131 42.04634922 39.4809524 37.21736696 35.20529102

33.40501255 31.78476192 30.31882087 28.9861472 27.76935819 26.65396827

25.62780953 24.68058609 23.80352735 22.98911566]

The first three terms, by observation, aren't linear.

Thank you so much for your hard work though!

[David Goodmanson]

Hi Aden,

C and D are constants in k, and depend on n only because they are a large-n limit. The expression

a_k = 1 + D/k    k>n/2

is definitely exact for the right value of D (which probably depends on n but has a large-n limit). Of course the method used in the answer did not get D quite correctly.

For the almost-linear region, If there were no n-dependent upper limit restrictions on the sum so that

a_k = 1 + (1/2k) Sum{j = 1,2k)} a_j

then it's not hard to show that

a_k = nC -2(k-1)

is an exact solution for all k. But it's not an exclusive solution and as I mentioned in the answer, the upper limit conditions set by n have an effect propagating down to smaller k.

I went back and modified the answer to make it more clear that the linear region is not exactly linear. Below is a plot of diff(a) and it would be great if it had a constant value of -2 for 1<k<n/2 but that's not the case.

The value of a_k depends on a sum up to 2k, and conversely any funky behavior at 2k could have an effect on a_k. In the plot I think the peak at n/4 might be due to the big discontinuity at n/2. If you run the plot an zoom in, you can better see another peak at n/8 (actually 130 not 125) and another at 72, not really n/16.

n = 1000;
m = zeros(n,n);
for k = 1:n
    q = min(n,2*k);
    m(k,1:q) = 1;
end
k = (1:n)';
m = diag(2*k)-m;
rhs = 2*k;
a = m\rhs;
% check the original set of equatons
chk = [];
for kk = 1:n
    q = min(n,2*kk);
    chk = [chk; [2*kk*a(kk) - (sum(a(1:q))+2*kk)]];
end
chkval = max(abs(chk))        % should be small
an = a/n;
C = (3/2-log(2))/(1-log(2))
D = (n/2)*(n*(C-1)+1)
k1 = (1:n/2)'; 
k2 = (n/2+1:n)';
b1 = n*C-2*(k1-1);
b2 = 1 + D./k2;
figure(1); grid on
plot(k,a/n,k1,b1/n,k2,b2/n)  % shown before in the original answer
ylim([0 3])
figure(2); grid on
plot(diff(a))

[David Goodmanson]

HI Aden, I modified the answer and also the comment above to make it clearer that the almost-linear part of the plot for 1<k<n/2 is not really linear.

[Aden]

Hi, I see. I understand your explanation and am convinced by it. Thank you for your detailed explanation!

However, one thing holding me back is the experimental values. Somehow, it really looks like it converges to the experimental 2.517782 instead of the theoretical 2.6294. Here is the database: https://docs.google.com/spreadsheets/d/1K-MNti6u8QtdoZlhCqj-AKb8KngqLhrSU9VpvI1Cmr0/edit?usp=sharing

Instead of directly finding the ratio , I figured that a better way was to find the difference between the terms (as in, ( for n equations) - ( for n-1 equations)). This way, it would converge faster.

[David Goodmanson]

Hi Aden,

The model value of 2.63 is only an approximate value, which is not too bad since it is within 4.5% of the real result, which is probably very close to 2.5177.... The 'theoretical' model value It was never intended to be taken as the actual answer. That's because model presumes the lower half of the plot to be linear, wheras the real calculation shows that it is close to, but actually not, linear. So there will be correction terms to the model, but I don't know how to obtain them. Your difference idea definitely seems worth trying.

Answer 2

[John D'Errico]

(Edited to fix my error in extrapolating your equation system. I had not seen at first you were essentially adding two terms to each equation for increasing k.)

For ANY value of n, you have n-1 unknowns, and n-1 equations. But all you want to do is to compute a_1, disregarding the other unknowns. And that would seem reasonable. Well, maybe possible.

I can create the problem in matrix form, as such:

function [A,b] = Ab(N)
  [K,J] = ndgrid(1:N-1);
  A = diag(2*(1:(N-1))) - (2*K >= J);
  b = 2*(1:(N-1))';
end

As you can see, I combined all terms into one matrix A.

[A,b] = Ab(4)

A = 3×3

1 -1 0 -1 3 -1 -1 -1 5

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

b = 3×1

2 4 6

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The complete solution for n==4 is then

format rat
A\b

ans = 
      27/4     
      19/4     
       7/2     

First, we should notice that A is always full rank. Just by looking at A in general, that would seem clear, but a proof should not be impossible.

[A,b] = Ab(5)

A = 
       1             -1              0              0       
      -1              3             -1             -1       
      -1             -1              5             -1       
      -1             -1             -1              7       

b = 
       2       
       4       
       6       
       8       

A\b

ans = 
      51/5     
      41/5     
      29/5     
      23/5     

[A,b] = Ab(6);
A\b

ans = 
     911/76    
     759/76    
     149/19    
     233/38    
      97/19    

And that clearly works, at least for specific values of n to yield the same result you have.

One idea might be to look at the LU factors of A. Note they are quite well behaved as n grows. For example:

[A,b] = Ab(7);[L,U] = lu(A)

L = 
       1              0              0              0              0              0       
      -1              1              0              0              0              0       
      -1             -1              1              0              0              0       
      -1             -1             -1/2            1              0              0       
      -1             -1             -1/2           -3/5            1              0       
      -1             -1             -1/2           -3/5           -6/19           1       

U = 
       1             -1              0              0              0              0       
       0              2             -1             -1              0              0       
       0              0              4             -2             -1             -1       
       0              0              0              5             -3/2           -3/2     
       0              0              0              0             38/5          -12/5     
       0              0              0              0              0            168/19    

[A,b] = Ab(8);[L,U] = lu(A)

L = 
       1              0              0              0              0              0              0       
      -1              1              0              0              0              0              0       
      -1             -1              1              0              0              0              0       
      -1             -1             -1/2            1              0              0              0       
      -1             -1             -1/2           -3/5            1              0              0       
      -1             -1             -1/2           -3/5           -6/19           1              0       
      -1             -1             -1/2           -3/5           -6/19          -5/14           1       

U = 
       1             -1              0              0              0              0              0       
       0              2             -1             -1              0              0              0       
       0              0              4             -2             -1             -1              0       
       0              0              0              5             -3/2           -3/2           -1       
       0              0              0              0             38/5          -12/5           -8/5     
       0              0              0              0              0            168/19         -40/19    
       0              0              0              0              0              0             78/7     

Noting that as n grows, it is only the last column of L and U that seem to change, I might not be surprised if one could write those factors down in some simple form. I don't see how at the moment, but it might be worth looking at. And you can recover the solution as

U\(L\b)

ans = 
    7939/468   
    7003/468   
     235/18    
     424/39    
     347/39    
     887/117   
     259/39    

A nice thing about LU factors is, if you are careful, you can get the last element of the solution, knowing only L, and the last diagonal element of U. And if we flip the matrix A from left to right, which is equivalent to resequencing the unknowns in reverse order....

[L,U] = lu(fliplr(A));
C = L\b;
C(end)/U(end,end)

ans = 
    7939/468   

I hope you see where I am going. If we could write down the LU factors, in a simple way as a function of N, for the flipped (left to right) matrix A, then solving for the FIRST element of a becomes simple. All you need to know is a way to compute the factro L, as well as the last diagonal element of U. A very nice thing about this is it might also give you a way to prove the long term (linear with n) behavior of the first element of a.

[Aden]

Hi! I'm not sure that your code here matches the sum in the question.

"A_n = @(n) diag(2*(1:(n-1))) - tril(ones(n-1),1);

b_n = @(n) 2*(1:(n-1))';"

However, is this matrix correct for A?

For n=5, the matrix produced should be

1 -1 0 0

-1 3 -1 -1

-1 -1 5 -1

-1 -1 -1 7

For row k < n/2, the number of "-1"s above the main diagonal should be increasing by 2 every row you go down.

(For n=5, a1 should be 51/5. For n=6, a1 should be 911/76.)

As such, your answer limit 7.2668730... is not the same as mine of 2.51778...

Thank you for your answer regardless though!

[John D'Errico]

No problem. I was in a hurry, so I'll check things over. I hope perhaps my re-written answer may give you a direction to follow.

Answer 3

[Torsten]

Here is my code, that provies a range of answers from a to b, inclusive.

Besides that you shouldn't give the same names to the lower bound of your range and the symbolic array a, what is your question ?

I used the below code, and it seems to give similar results as yours. The behaviour of a1 with n looks quite linear - that's why I made a linear fit at the end.

nmax = 25;
a1 = sym('a1',[nmax,1]);
for n = 1:nmax
  a = sym('a',[n 1]);
  for k = 1:n
    ub = min(n-1,2*k);
    eqn(k) = 2*k*a(k)==sum(a(1:ub))+2*k;
  end
  [A,b] = equationsToMatrix(eqn,a);
  sol = A\b;
  a1(n) = sol(1);
end
a1 = double(a1);
M = [ones(nmax,1),(1:nmax).'];
rhs = a1;
sol = M\rhs;
hold on
plot(1:nmax,a1)
plot(1:nmax,sol(1)+sol(2)*(1:nmax))
hold off
grid on

[Aden]

Thank you!

Edit: In fact a_1/n should roughly converge to the number 2.51778 as n increases, (but I cannot prove it).

My question is, to find (if possible) a closed form expression of a_1 for each n. (In other words, a general solution for each n).

I get that this is slightly off-topic in relation to coding, but perhaps that this linear equation may be solved with row operations in some ingenious way with Cramer's rule, but I couldn't make any progress on that, so I resulted to coding.

[Torsten]

I have no idea. But let's wait for @John D'Errico - he promised:

I'll post an answer now, first, in terms of how I would solve it for specific n, and then look to see how we might generalize it for any n.

My question is, to find (if possible) a closed form expression of a_1 for each n. (In other words, a general solution for each n).

I don't fully understand what you mean here. I thought of a_1 as a function of n, like a_1 = 2.51778*n + sin(n^2). Or do you mean something else here ?

[John D'Errico]

Still writing. It will take an extra hour though, as my wife is asking for a boat ride. ;-) But yes, my hope is to see if a solution exists as a function of n, and possibly as a limit.

[Torsten]

@Aden

May I ask about the context in which this problem appeared ?

Is this a special field of investigation ?

Try asking the question in a pure maths forum:

https://math.stackexchange.com/