How to solve a linear ill conditioned system of equations
36 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
KostasK
el 7 de Nov. de 2021
Comentada: Matt J
el 8 de Nov. de 2021
Hi all, I have a typical linear problem to solve
, where A is a tridiagonal, illconditioned matrix. Using the solution of said system, I calculate a vector R, such that R = kv.*diff(x). I have what I know to be the correct value of R, and in my subsequent attempts I perform this calculation involving the solution that I obtain of the linear system, and compare to this value of R that I have.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/793299/image.png)
I use two ways to work around the ill-conditioning:
- I use the typical L-U decomposition, however I still get a warning and the results are not that accurate as you can see.
- I perform the same process as 1, however I condition U by adding the identity matrix which appears to solve the ill conditioning problem. Moreover, I notice that the result R2 is away by the constant R2(1) from the known answer R, so by correcting for that this approach appears to give me the right answer.
clear ; clc
A = [663100000,-663100000,0,0,0,0,0,0,0,0,0,0;-663100000,1407100000.00000,-744000000,0,0,0,0,0,0,0,0,0;0,-744000000,1488000000.00000,-744000000,0,0,0,0,0,0,0,0;0,0,-744000000,1488000000.00000,-744000000,0,0,0,0,0,0,0;0,0,0,-744000000,1488000000.00000,-744000000,0,0,0,0,0,0;0,0,0,0,-744000000,1488000000.00000,-744000000,0,0,0,0,0;0,0,0,0,0,-744000000,1865000000.00000,-1121000000.00000,0,0,0,0;0,0,0,0,0,0,-1121000000.00000,3086000000.00000,-1965000000.00000,0,0,0;0,0,0,0,0,0,0,-1965000000.00000,51965000000.0000,-50000000000.0000,0,0;0,0,0,0,0,0,0,0,-50000000000.0000,50022150000.0000,-22150000,0;0,0,0,0,0,0,0,0,0,-22150000,111860000,-89710000;0,0,0,0,0,0,0,0,0,0,-89710000,89710000] ;
b = [0;97346.0079330557;97020.9494874157;97246.0491601032;97462.7637751785;97128.3959227217;97326.7056988805;0;0;0;0;-460875.388953404] ;
R = [0;97346.0079330557;194366.957420471;291613.006580575;389075.770355753;486204.166278475;583530.871977355;583530.871977355;583530.871977355;583530.871977355;583530.871977355] ;
kv = diag(A, -1) ;
% L U Decomposition
[L, U] = lu(A) ;
x1 = U\(L\b) ;
R1 = kv.*diff(x1) ;
% L U Decomposition with Preconditioning
x2 = (U + eye(length(A)))\(L\b) ;
R2 = kv.*diff(x2) ;
% Compare results. R is the correct solution of the problem
[R R1 R2-R2(1)]
Based on the above I have two questions:
- I cannot explain why my second approach works, and specifically why I have to perform the calculation R2-R2(1) in order to get the right result. Is this a coincidence for this case, or can I safely generalise this for the other similar ill-conditioned matrices like A that I have so that I can solve similar problems correctly?
- In case that the above is not a reliable solution, can you suggest of a method to handle this problem? (I have looked into the Tikhonov regularisation in regtools however I cannot get the function to work as it keeps returning NaN.)
Thanks for your help in advance.
2 comentarios
Matt J
el 7 de Nov. de 2021
I have a typical linear problem to solve , where A is a tridiagonal, illconditioned matrix. I have the correct answer to be R = kv.*diff(b)
That clearly isn't the case in your example. R does not solve the original equations:
A = [663100000,-663100000,0,0,0,0,0,0,0,0,0,0;-663100000,1407100000.00000,-744000000,0,0,0,0,0,0,0,0,0;0,-744000000,1488000000.00000,-744000000,0,0,0,0,0,0,0,0;0,0,-744000000,1488000000.00000,-744000000,0,0,0,0,0,0,0;0,0,0,-744000000,1488000000.00000,-744000000,0,0,0,0,0,0;0,0,0,0,-744000000,1488000000.00000,-744000000,0,0,0,0,0;0,0,0,0,0,-744000000,1865000000.00000,-1121000000.00000,0,0,0,0;0,0,0,0,0,0,-1121000000.00000,3086000000.00000,-1965000000.00000,0,0,0;0,0,0,0,0,0,0,-1965000000.00000,51965000000.0000,-50000000000.0000,0,0;0,0,0,0,0,0,0,0,-50000000000.0000,50022150000.0000,-22150000,0;0,0,0,0,0,0,0,0,0,-22150000,111860000,-89710000;0,0,0,0,0,0,0,0,0,0,-89710000,89710000] ;
b = [0;97346.0079330557;97020.9494874157;97246.0491601032;97462.7637751785;97128.3959227217;97326.7056988805;0;0;0;0;-460875.388953404] ;
R = [0;97346.0079330557;194366.957420471;291613.006580575;389075.770355753;486204.166278475;583530.871977355;583530.871977355;583530.871977355;583530.871977355;583530.871977355] ;
A*R-b
Respuesta aceptada
Matt J
el 7 de Nov. de 2021
Editada: Matt J
el 8 de Nov. de 2021
I don't think the equations should require regularization. The only vector in the nullspace of A is ones(12,1) , but R is insensitive to this null vector because kv.*(diff(x+c*ones(12,1))= kv.*diff(x).
The discrepancies you are getting in version R1 seem to be either because your expected values for R are incorrect or else there is something wrong in your A,b data, probably the final row. I suspect this because of the results I get below from lsqlin. Here I am solving for the least squares solution to A*x=b subject to the constraint R=kv.*diff(x). As you can see, the final equation in A*x=b is not solved very well compared to the other equations.
A = [663100000,-663100000,0,0,0,0,0,0,0,0,0,0;-663100000,1407100000.00000,-744000000,0,0,0,0,0,0,0,0,0;0,-744000000,1488000000.00000,-744000000,0,0,0,0,0,0,0,0;0,0,-744000000,1488000000.00000,-744000000,0,0,0,0,0,0,0;0,0,0,-744000000,1488000000.00000,-744000000,0,0,0,0,0,0;0,0,0,0,-744000000,1488000000.00000,-744000000,0,0,0,0,0;0,0,0,0,0,-744000000,1865000000.00000,-1121000000.00000,0,0,0,0;0,0,0,0,0,0,-1121000000.00000,3086000000.00000,-1965000000.00000,0,0,0;0,0,0,0,0,0,0,-1965000000.00000,51965000000.0000,-50000000000.0000,0,0;0,0,0,0,0,0,0,0,-50000000000.0000,50022150000.0000,-22150000,0;0,0,0,0,0,0,0,0,0,-22150000,111860000,-89710000;0,0,0,0,0,0,0,0,0,0,-89710000,89710000] ;
b = [0;97346.0079330557;97020.9494874157;97246.0491601032;97462.7637751785;97128.3959227217;97326.7056988805;0;0;0;0;-460875.388953404] ;
R = [0;97346.0079330557;194366.957420471;291613.006580575;389075.770355753;486204.166278475;583530.871977355;583530.871977355;583530.871977355;583530.871977355;583530.871977355] ;
s=1e6;
A=A/s; b=b/s; R=R/s;
kv = diag(A, -1) ;
[m,n]=size(A);
E=kv(:).*diff(eye(n));
opts=optimoptions('lsqlin','OptimalityTolerance',1e-20,'ConstraintTolerance',1e-20);
[x,fval]=lsqlin(A,b,[],[],E,R ,[],[],[],opts); fval
equationError=A*x-b;
equationError([1:4,10:12])
2 comentarios
Matt J
el 8 de Nov. de 2021
I don't think you need to. I think you can just do the calculation as,
R=kv.*diff(pinv(A)*b)
Más respuestas (0)
Ver también
Categorías
Más información sobre Linear Algebra en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!