Generalised solution to a n-order system of equations
97 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Aden
am 11 Sep. 2024
Bearbeitet: David Goodmanson
am 15 Sep. 2024 um 20:41
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.
7 Kommentare
John D'Errico
am 12 Sep. 2024 um 14:22
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.
Akzeptierte Antwort
David Goodmanson
am 13 Sep. 2024 um 20:34
Bearbeitet: David Goodmanson
am 14 Sep. 2024 um 20:22
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%.
5 Kommentare
David Goodmanson
am 15 Sep. 2024 um 20:38
Bearbeitet: David Goodmanson
am 15 Sep. 2024 um 20:41
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.
Weitere Antworten (2)
John D'Errico
am 12 Sep. 2024 um 14:52
Bearbeitet: John D'Errico
am 12 Sep. 2024 um 18:35
(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)
The complete solution for n==4 is then
format rat
A\b
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\b
[A,b] = Ab(6);
A\b
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)
[A,b] = Ab(8);[L,U] = lu(A)
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)
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)
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.
2 Kommentare
John D'Errico
am 12 Sep. 2024 um 16:43
Bearbeitet: John D'Errico
am 12 Sep. 2024 um 18:31
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.
Torsten
am 12 Sep. 2024 um 14:20
Bearbeitet: Torsten
am 12 Sep. 2024 um 14:21
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
4 Kommentare
John D'Errico
am 12 Sep. 2024 um 14:55
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.
Siehe auch
Kategorien
Mehr zu Polynomials finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!