Matrices. Sistema de ecuaciones lineales. Valores y vectores propios.
Determinante y matriz inversa de una matriz cuadrada
Calcular el determinate y la matriz inversa de.
Comprueba que A·A-1=I (matriz identidad)
Solución
public class Vector {
public int n; //dimensión
double[] x;
public Vector(int n) {
this.n=n;
x=new double[n];
for(int i=0; i<n; i++){
x[i]=0.0;
}
}
public Vector(double[] x) {
this.x=x;
n=x.length;
}
public String toString(){
String texto=" ";
for(int i=0; i<n; i++){
texto+="\t "+(double)Math.round(1000*x[i])/1000;
}
texto+="\n";
return texto;
}
public static Vector producto(Matriz a,Vector v){
int n=v.n;
Vector b=new Vector(n);
for(int i=0; i<n; i++){
for(int k=0; k<n; k++){
b.x[i]+=a.x[i][k]*v.x[k];
}
}
return b;
}
}
public class Matriz implements Cloneable{
public int n; //dimensión
public double[][] x;
public Matriz(int n) {
this.n=n;
x=new double[n][n];
for(int i=0; i<n; i++){
for(int j=0; j<n; j++){
x[i][j]=0.0;
}
}
}
public Matriz(double[][] x) {
this.x=x;
n=x.length;
}
public String toString(){
String texto="\n";
for(int i=0; i<n; i++){
for(int j=0; j<n; j++){
texto+="\t "+(double)Math.round(1000*x[i][j])/1000;
}
texto+="\n";
}
texto+="\n";
return texto;
}
public Object clone(){
Matriz obj=null;
try{
//llama a clone de la clase base Object
obj=(Matriz)super.clone();
}catch(CloneNotSupportedException ex){
System.out.println(" no se puede duplicar");
}
//copia la matriz bidimensional
obj.x=(double[][])obj.x.clone();
for(int i=0; i<obj.x.length; i++){
obj.x[i]=(double[])obj.x[i].clone();
}
return obj;
}
public static Matriz producto(Matriz a, Matriz b){
Matriz resultado=new Matriz(a.n);
for(int i=0; i<a.n; i++){
for(int j=0; j<a.n; j++){
for(int k=0; k<a.n; k++){
resultado.x[i][j]+=a.x[i][k]*b.x[k][j];
}
}
}
return resultado;
}
public double determinante(){
Matriz a=(Matriz)clone();
for(int k=0; k<n-1; k++){
for(int i=k+1; i<n; i++){
for(int j=k+1; j<n; j++){
a.x[i][j]-=a.x[i][k]*a.x[k][j]/a.x[k][k];
}
}
}
double deter=1.0;
for(int i=0; i<n; i++){
deter*=a.x[i][i];
}
return deter;
}
public static Matriz inversa(Matriz d){
int n=d.n;
Matriz a=(Matriz)d.clone();
Matriz b=new Matriz(n);
Matriz c=new Matriz(n);
//matriz unidad
for(int i=0; i<n; i++){
b.x[i][i]=1.0;
}
//transformación de la matriz y de los términos independientes
for(int k=0; k<n-1; k++){
for(int i=k+1; i<n; i++){
//términos independientes
for(int s=0; s<n; s++){
b.x[i][s]-=a.x[i][k]*b.x[k][s]/a.x[k][k];
}
//elementos de la matriz
for(int j=k+1; j<n; j++){
a.x[i][j]-=a.x[i][k]*a.x[k][j]/a.x[k][k];
}
}
}
//cálculo de las incógnitas, elementos de la matriz inversa
for(int s=0; s<n; s++){
c.x[n-1][s]=b.x[n-1][s]/a.x[n-1][n-1];
for(int i=n-2; i>=0; i--){
c.x[i][s]=b.x[i][s]/a.x[i][i];
for(int k=n-1; k>i; k--){
c.x[i][s]-=a.x[i][k]*c.x[k][s]/a.x[i][i];
}
}
}
return c;
}
}
public class MyClass1 {
public static void main(String[] args){
double[][] a1={{3,1,-1,2,1},{-2,3,1,4,3},{1,4,2,3,1},{5,-2,-3,5,-1},{-1,1,2,3,2}};
Matriz a=new Matriz(a1);
System.out.println("Determinante "+a.determinante());
System.out.println("Inversa: "+Matriz.inversa(a));
System.out.println("Comprobación:"+Matriz.producto(a, Matriz.inversa(a)));
}
}
Sistema de ecuaciones lineales
Resolver el sistema de eucaciones lineales
Solución
//Poner aquí las definiciones de las clases Vector y Matriz (ejercicio 1)
public class MyClass1 {
public static void main(String[] args){
double[][] m={{Math.PI,-Math.E,Math.sqrt(2),-Math.sqrt(3)},
{Math.PI*Math.PI, Math.E,-Math.E*Math.E,3.0/7},
{Math.sqrt(5),-Math.sqrt(6),1,-Math.sqrt(2)},
{Math.PI*Math.PI*Math.PI,Math.E*Math.E,-Math.sqrt(7),1.0/9}};
Matriz coef=new Matriz(m);
double[] n={Math.sqrt(11),0,Math.PI,Math.sqrt(2)};
Vector ter=new Vector (n);
Vector solucion=(Vector.producto(Matriz.inversa(coef), ter));
System.out.println("solución :"+solucion);
}
}
Resolver el sistema de ecuaciones lineales
Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P. Numerical Recipes in C, Second edition. Chapter 2. Solution of Linear Algebraic Equations. 2.3 LU Decomposition and Its Applications. Cambridge University Press. Código en C adaptado por el autor al lenguaje Java
Solución
public class LU {
private final static double TINY=1.0e-20;
static private void ludcmp(double[][] a, int[] indx) {
int n=a[0].length;
int i,imax=0,j,k;
double big,dum,sum,temp;
double[] vv=new double[n];
for (i=0;i<n;i++) {
big=0.0;
for (j=0;j<n;j++)
if ((temp=Math.abs(a[i][j])) > big) big=temp;
if (big == 0.0) System.out.println("Singular matrix in routine ludcmp");
vv[i]=1.0/big;
}
for (j=0;j<n;j++) {
for (i=0;i<j;i++) {
sum=a[i][j];
for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<n;i++) {
sum=a[i][j];
for (k=0;k<j;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (dum=vv[i]*Math.abs(sum)) >= big) {
big=dum;
imax=i;
}
}
if (j != imax) {
for (k=0;k<n;k++) {
dum=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
}
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j][j] == 0.0) a[j][j]=TINY;
if (j != n-1) {
dum=1.0/(a[j][j]);
for (i=j+1;i<n;i++) a[i][j] *= dum;
}
}
}
static private void lubksb(double[][] a, double b[], int[] indx){
int n=a[0].length;
int i,ii=0,ip,j;
double sum;
for (i=1;i<=n;i++) {
ip=indx[i-1];
sum=b[ip];
b[ip]=b[i-1];
if (ii>0)
for (j=ii;j<=i-1;j++) sum -= a[i-1][j-1]*b[j-1];
else if (sum!=0) ii=i;
b[i-1]=sum;
}
for (i=n-1;i>=0;i--) {
sum=b[i];
for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
b[i]=sum/a[i][i];
}
}
static public void resolver(double[][]a, double[] b){
int[] index=new int[b.length];
LU.ludcmp(a, index);
LU.lubksb(a, b, index);
}
}
public class MyClass1 {
public static void main(String[] args) {
double[][] a={{3, 1, -1, 2}, {-5, 1, 3, -4}, {2, 0, 1, -1}, {1, -5, -3, -3}};
double[] b={6, -12, 1, 3};
LU.resolver(a, b);
//solución
for(int i=0; i<b.length; i++){
System.out.print(b[i]+" ");
}
}
}
Valores propios y vectores propios
Calcular los valores propios de la matriz simétrica
Solución
public class Vector {
public int n;
double[] x;
public Vector(int n) {
this.n=n;
x=new double[n];
for(int i=0; i<n; i++){
x[i]=0.0;
}
}
public Vector(double[] x) {
this.x=x;
n=x.length;
}
public String toString(){
String texto=" ";
for(int i=0; i<n; i++){
texto+="\t "+(double)Math.round(1000*x[i])/1000;
}
texto+="\n";
return texto;
}
public static Vector producto(Matriz a,Vector v){
int n=v.n;
Vector b=new Vector(n);
for(int i=0; i<n; i++){
for(int k=0; k<n; k++){
b.x[i]+=a.x[i][k]*v.x[k];
}
}
return b;
}
}
public class Matriz implements Cloneable{
public int n;
public double[][] x;
public Matriz(int n) {
this.n=n;
x=new double[n][n];
for(int i=0; i<n; i++){
for(int j=0; j<n; j++){
x[i][j]=0.0;
}
}
}
public Matriz(double[][] x) {
this.x=x;
n=x.length;
}
public String toString(){
String texto="\n";
for(int i=0; i<n; i++){
for(int j=0; j<n; j++){
texto+="\t "+(double)Math.round(1000*x[i][j])/1000;
}
texto+="\n";
}
texto+="\n";
return texto;
}
public Object clone(){
Matriz obj=null;
try{
//llama a clone de la clase base Object
obj=(Matriz)super.clone();
}catch(CloneNotSupportedException ex){
System.out.println(" no se puede duplicar");
}
//copia la matriz bidimensional
obj.x=(double[][])obj.x.clone();
for(int i=0; i<obj.x.length; i++){
obj.x[i]=(double[])obj.x[i].clone();
}
return obj;
}
public static Matriz producto(Matriz a, Matriz b){
Matriz resultado=new Matriz(a.n);
for(int i=0; i<a.n; i++){
for(int j=0; j<a.n; j++){
for(int k=0; k<a.n; k++){
resultado.x[i][j]+=a.x[i][k]*b.x[k][j];
}
}
}
return resultado;
}
public static Matriz inversa(Matriz d){
int n=d.n;
Matriz a=(Matriz)d.clone();
Matriz b=new Matriz(n);
Matriz c=new Matriz(n);
//matriz unidad
for(int i=0; i<n; i++){
b.x[i][i]=1.0;
}
//transformación de la matriz y de los términos independientes
for(int k=0; k<n-1; k++){
for(int i=k+1; i<n; i++){
//términos independientes
for(int s=0; s<n; s++){
b.x[i][s]-=a.x[i][k]*b.x[k][s]/a.x[k][k];
}
//elementos de la matriz
for(int j=k+1; j<n; j++){
a.x[i][j]-=a.x[i][k]*a.x[k][j]/a.x[k][k];
}
}
}
//cálculo de las incógnitas, elementos de la matriz inversa
for(int s=0; s<n; s++){
c.x[n-1][s]=b.x[n-1][s]/a.x[n-1][n-1];
for(int i=n-2; i>=0; i--){
c.x[i][s]=b.x[i][s]/a.x[i][i];
for(int k=n-1; k>i; k--){
c.x[i][s]-=a.x[i][k]*c.x[k][s]/a.x[i][i];
}
}
}
return c;
}
//matriz traspuesta
static Matriz traspuesta(Matriz a){
int n=a.n;
Matriz d=new Matriz(a.n);
for(int i=0; i<n; i++){
for(int j=0; j<n; j++){
d.x[i][j]=a.x[j][i];
}
}
return d;
}
public double[] polCaracteristico(){
Matriz pot=new Matriz(n);
//matriz unidad
for(int i=0; i<n; i++){
pot.x[i][i]=1.0;
}
double[] p=new double[n+1];
double[] s=new double[n+1];
//potencias de la matriz y traza
for(int i=1; i<=n; i++){
pot=Matriz.producto(pot, this);
s[i]=pot.traza();
}
//coeficientes del polinomio característico
p[0]=1.0;
p[1]=-s[1];
for(int i=2; i<=n; i++){
p[i]=-s[i]/i;
for(int j=1; j<i; j++){
p[i]-=s[i-j]*p[j]/i;
}
}
return p;
}
public double traza(){
double t=0;
for(int i=0; i<n; i++) t+=x[i][i];
return t;
}
}
//Raíces del polinomio por el procedimiento de Graeffe
//Vectores propios
//Definición de las clase Graeffe y Complejo
public class MyClass1 {
public static void main(String[] args) {
double[][] p1={{1,2,3,4},{2,1,2,3},{3,2,1,2},{4,3,2,1}};
Matriz p=new Matriz(p1);
double pol[]=p.polCaracteristico();
System.out.println("Polinomio característico");
for(int i=0; i<pol.length; i++){
System.out.print((double)Math.round(pol[i]*1000)/1000+" , ");
}
System.out.println("");
//raíces del polinomio
Graeffe g=new Graeffe(pol);
g.mostrarRaices();
double[] lambda=g.raicesReales; //valores propios
//vectores propios
double[][] p2=new double[p1.length-1][p1.length-1];
double[] ind=new double[p2.length];
Vector x;
Matriz m;
for(int j=0; j<ind.length; j++){
ind[j]=-p1[0][j+1];
}
Vector v=new Vector(ind);
for(int k=0; k<lambda.length; k++){
for(int i=0; i<p2.length; i++){
for(int j=0; j<p2.length; j++){
p2[i][j]=p1[i+1][j+1];
if(i==j) p2[i][j]=p1[i+1][j+1]-lambda[k];
}
}
m=new Matriz(p2);
System.out.println("Valor propio "+ lambda[k]);
System.out.println("Vectores propios ");
x=Vector.producto(Matriz.inversa(m), v);
System.out.println("1");
System.out.println(x);
}
}
}
Calcular los valores y vectores propios de la matriz simétrica por el procedimiento de Jacobi
Solución
public class Matriz implements Cloneable{
public int n; //dimensión
private double[][] x;
public Matriz(int n) {
this.n=n;
x=new double[n][n];
for(int i=0; i<n; i++){
for(int j=0; j<n; j++){
x[i][j]=0.0;
}
}
}
public Matriz(double[][] x) {
this.x=x;
n=x.length;
}
public Object clone(){
Matriz obj=null;
try{
obj=(Matriz)super.clone();
}catch(CloneNotSupportedException ex){
System.out.println(" no se puede duplicar");
}
//este es el código para clonar la matriz bidimensional
obj.x=(double[][])obj.x.clone();
for(int i=0; i<obj.x.length; i++){
obj.x[i]=(double[])obj.x[i].clone();
}
return obj;
}
//producto de dos matrices
static Matriz producto(Matriz a, Matriz b){
Matriz resultado=new Matriz(a.n);
for(int i=0; i<a.n; i++){
for(int j=0; j<a.n; j++){
for(int k=0; k<a.n; k++){
resultado.x[i][j]+=a.x[i][k]*b.x[k][j];
}
}
}
return resultado;
}
//producto de una matriz por un escalar
static Matriz producto(Matriz a, double d){
Matriz resultado=new Matriz(a.n);
for(int i=0; i<a.n; i++){
for(int j=0; j<a.n; j++){
resultado.x[i][j]=a.x[i][j]*d;
}
}
return resultado;
}
//producto de un escalar por una matriz
static Matriz producto(double d, Matriz a){
Matriz resultado=new Matriz(a.n);
for(int i=0; i<a.n; i++){
for(int j=0; j<a.n; j++){
resultado.x[i][j]=a.x[i][j]*d;
}
}
return resultado;
}
//matriz traspuesta
static Matriz traspuesta(Matriz a){
int n=a.n;
Matriz d=new Matriz(a.n);
for(int i=0; i<n; i++){
for(int j=0; j<n; j++){
d.x[i][j]=a.x[j][i];
}
}
return d;
}
public Matriz valoresPropios(double[] valores, int maxIter)throws ValoresExcepcion{
final double CERO=1e-8;
double maximo, tolerancia, sumsq;
double x, y, z, c, s;
int contador=0;
int i, j, k, l;
Matriz a=(Matriz)clone(); //matriz copia
Matriz p=new Matriz(n);
Matriz q=new Matriz(n);
//matriz unidad
for(i=0; i<n; i++){
q.x[i][i]=1.0;
}
do{
k=0; l=1;
maximo=Math.abs(a.x[k][1]);
for(i=0; i<n-1; i++){
for(j=i+1; j<n; j++){
if(Math.abs(a.x[i][j])>maximo){
k=i; l=j;
maximo=Math.abs(a.x[i][j]);
}
}
}
sumsq=0.0;
for(i=0; i<n; i++){
sumsq+=a.x[i][i]*a.x[i][i];
}
tolerancia=0.0001*Math.sqrt(sumsq)/n;
if(maximo<tolerancia) break;
//calcula la matriz ortogonal de p
//inicialmente es la matriz unidad
for(i=0; i<n; i++){
for(j=0; j<n; j++){
p.x[i][j]=0.0;
}
}
for(i=0; i<n; i++){
p.x[i][i]=1.0;
}
y=a.x[k][k]-a.x[l][l];
if(Math.abs(y)<CERO){
c=s=Math.sin(Math.PI/4);
}else{
x=2*a.x[k][l];
z=Math.sqrt(x*x+y*y);
c=Math.sqrt((z+y)/(2*z));
s=signo(x/y)*Math.sqrt((z-y)/(2*z));
}
p.x[k][k]=c;
p.x[l][l]=c;
p.x[k][l]=s;
p.x[l][k]=-s;
a=Matriz.producto(p, Matriz.producto(a, Matriz.traspuesta(p)));
q=Matriz.producto(q, Matriz.traspuesta(p));
contador++;
}while(contador<maxIter);
if(contador==maxIter){
throw new ValoresExcepcion("No se han podido calcular los valores propios");
}
//valores propios
for(i=0; i<n; i++){
valores[i]=(double)Math.round(a.x[i][i]*1000)/1000;
}
//vectores propios
return q;
}
private int signo(double x){
return (x>0 ? 1 : -1);
}
public String toString(){
String texto="\n";
for(int i=0; i<n; i++){
for(int j=0; j<n; j++){
texto+="\t "+(double)Math.round(1000*x[i][j])/1000;
}
texto+="\n";
}
texto+="\n";
return texto;
}
}
class ValoresExcepcion extends Exception {
public ValoresExcepcion() {
super();
}
public ValoresExcepcion(String s) {
super(s);
}
}public class MyClass1 {
public static void main(String[] args) {
System.out.println("Valores y vectores propios ");
double[][] p1={{1,2,3,4},{2,1,2,3},{3,2,1,2},{4,3,2,1}};
Matriz matriz=new Matriz(p1);
Matriz vectores=new Matriz(matriz.n);
double[] valores=new double[matriz.n];
try{
vectores=matriz.valoresPropios(valores, 20);
}catch(ValoresExcepcion ex){
System.out.println("Al calcular los valores propios
se ha producido una excepción\n "
+ex.getClass()+ " con el mensaje "+ ex.getMessage());
}
for(int i=0; i<matriz.n; i++){ //valores propios
System.out.print(valores[i]+" , ");
}
System.out.println("");
System.out.println(vectores); //vectores propios
}
}
Calcular los valores y vectores propios de la matriz simétrica
Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P. Numerical Recipes in C, Second edition. Chapter 11. Eigensystems. 11.1 Jacobi Transformations of a Symmetric Matrix. Cambridge University Press. Código en C adaptado por el autor al lenguaje Java
Solución
public class Jacobi {
//el vector d son los valores propios
//Las columnas de la matriz v son los vectores propios
//devuelve el número de iteracciones
static int calcula(double[][] a, double[] d, double[][] v) {
int n=a[0].length;
int j,iq,ip,i;
double tresh,theta,tau,t,sm,s,h,g,c;
double[] b=new double[n]; //b=vector(1,n);
double[] z=new double[n]; //z=vector(1,n);
for (ip=0;ip<n;ip++) { //Initialize to the identity matrix.
for (iq=0;iq<n;iq++) v[ip][iq]=0.0;
v[ip][ip]=1.0;
}
for (ip=0;ip<n;ip++) { //Initialize b and d to the diagonal of a.
b[ip]=d[ip]=a[ip][ip];
z[ip]=0.0;
//This vector will accumulate terms of the form tapq as in equation
}
int nrot=0;
for (i=1;i<=50;i++) {
sm=0.0;
for (ip=0;ip<=n-2;ip++) { //Sum o -diagonal elements.
for (iq=ip+1;iq<n;iq++)
sm += Math.abs(a[ip][iq]);
}
if (sm == 0.0) {
//The normal return, which relies on quadratic convergence to machine underflow.
return nrot;
}
if (i < 4)
tresh=0.2*sm/(n*n);// ...on the rst three sweeps.
else
tresh=0.0; //...thereafter.
for (ip=0;ip<n-1;ip++) {
for (iq=ip+1;iq<n;iq++) {
g=100.0*Math.abs(a[ip][iq]);
//After four sweeps, skip the rotation if the o -diagonal element is small.
if (i > 4 && (float)(Math.abs(d[ip])+g) == (float)Math.abs(d[ip])
&& (float)(Math.abs(d[iq])+g) == (float)Math.abs(d[iq]))
a[ip][iq]=0.0;
else if (Math.abs(a[ip][iq]) > tresh) {
h=d[iq]-d[ip];
if ((float)(Math.abs(h)+g) == (float)Math.abs(h))
t=(a[ip][iq])/h; // t = 1=(2 )
else {
theta=0.5*h/(a[ip][iq]); //Equation (11.1.10).
t=1.0/(Math.abs(theta)+Math.sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
c=1.0/Math.sqrt(1+t*t);
s=t*c;
tau=s/(1.0+c);
h=t*a[ip][iq];
z[ip] -= h;
z[iq] += h;
d[ip] -= h;
d[iq] += h;
a[ip][iq]=0.0;
for (j=0;j<=ip-1;j++) {
g=a[j][ip];
h=a[j][iq];
a[j][ip]=g-s*(h+g*tau);
a[j][iq]=h+s*(g-h*tau);
}
for (j=ip+1;j<=iq-1;j++) {
g=a[ip][j];
h=a[j][iq];
a[ip][j]=g-s*(h+g*tau);
a[j][iq]=h+s*(g-h*tau);
}
for (j=iq+1;j<n;j++) {
g=a[ip][j];
h=a[iq][j];
a[ip][j]=g-s*(h+g*tau);
a[iq][j]=h+s*(g-h*tau);
}
for (j=0;j<n;j++) {
g=v[j][ip];
h=v[j][iq];
v[j][ip]=g-s*(h+g*tau);
v[j][iq]=h+s*(g-h*tau);
}
++nrot;
}
}
}
for (ip=0;ip<n;ip++) {
b[ip] += z[ip];
d[ip]=b[ip]; //Update d with the sum of tapq,
z[ip]=0.0; //and reinitialize z.
}
}
return nrot;
}
}
public class MyClass1 {
public static void main(String[] args) {
System.out.println("Valores y vectores propios ");
double[][] m={{4,2,2},{2,5,1},{2,1,6}};
double[][] vectores=new double[m[0].length][m[0].length];
double[] valores=new double[m[0].length];
Jacobi.calcula(m, valores, vectores);
System.out.println("Valores propios");
for(int i=0; i<m[0].length; i++){
System.out.print(valores[i]+" , ");
}
System.out.println("");
System.out.println("Vectores propios");
for(int i=0; i<m[0].length; i++){
for(int j=0; j<m[0].length; j++){
System.out.print(vectores[i][j]+" , ");
}
System.out.println();
}
}
}