码迷,mamicode.com
首页 > 其他好文 > 详细

Commemoration of Linear Algebra

时间:2015-04-17 22:07:41      阅读:146      评论:0      收藏:0      [点我收藏+]

标签:

 

1. Gaussian Elimination

  1 import java.util.*;
  2 import java.text.*;
  3 
  4 class Linear {
  5     private int m, n;
  6     private double[][] mat;
  7     private DecimalFormat df;
  8     
  9     public Linear(int m,int n) {
 10         this.m = m; 
 11         this.n = n;
 12         mat = new double[m][n+1];
 13         df = new DecimalFormat("#.00");
 14     }
 15     public void getData(Scanner in) {
 16         for (int i=0;i<m;i++) {
 17             for (int j=0;j<=n;j++) {
 18                 mat[i][j] = in.nextDouble();
 19             }
 20         }
 21     }
 22     public void solve() {
 23         for (int j=0;j<m&&j<n;j++) {
 24             utilize(j);
 25             for (int i=j+1;i<m;i++) {
 26                 subtract(i,j);
 27             }
 28         }
 29         int r = rank();
 30         for (int i=r;i<m;i++) {
 31             if (!zero(mat[i][n])) {
 32                 System.out.println("NO SOLUTION");
 33                 return;
 34             }
 35         }
 36         for (int i=r-2;i>=0;i--) {
 37             for (int j=i+1;j<r;j++) {
 38                 subtract(i,j);
 39             }
 40         }
 41         printSol(r);
 42     }
 43     private void utilize(int j) {
 44         int t = m-1;
 45         for (int i=j;i<=t;i++) {
 46             if (zero(mat[i][j])) {
 47                 swapRow(i--,t--);
 48             } else {
 49                 for (int k=n;k>=j;k--) {
 50                     mat[i][k]/=mat[i][j];
 51                 }
 52             }
 53         }
 54     }
 55     private void swapRow(int i,int t) {
 56         for (int j=0;j<=n;j++) {
 57             double tmp = mat[i][j];
 58             mat[i][j] = mat[t][j];
 59             mat[t][j] = tmp;
 60         }
 61     }
 62     private int rank() {
 63         int r = 0;
 64         for (;r<m&&r<n;r++) {
 65             if (zero(mat[r][r])) {
 66                 break;
 67             }
 68         }
 69         return r;
 70     }
 71     private void subtract(int i,int j) {
 72         if (zero(mat[j][j])) {
 73             return;
 74         }
 75         for (int k=n;k>=j;k--) {
 76             mat[i][k]-=mat[i][j]*mat[j][k];
 77         }
 78     }
 79     private void printSol(int r) {
 80         System.out.print("X = \t");
 81         parSol(r);
 82         for (int j=r;j<n;j++) {
 83             genSol(r,j);
 84         }
 85     }
 86     private void parSol(int r) {
 87         System.out.print("(");
 88         if (r==0) {
 89             System.out.print(df.format(0));
 90             for (int i=1;i<n;i++) {
 91                 System.out.print(", "+df.format(0));
 92             }
 93         } else {
 94             System.out.print(df.format(mat[0][n]));
 95             for (int i=1;i<r;i++) {
 96                 System.out.print(", "+df.format(mat[i][n]));
 97             }
 98             for (int i=r;i<n;i++) {
 99                 System.out.print(", "+df.format(0));
100             }
101         }
102         System.out.println(")‘");
103     }
104     private void genSol(int r,int j) {
105         System.out.print(" + k"+(j-r+1)+" *\t(");
106         if (r==0) {
107             for (int i=0;i<j;i++) {
108                 System.out.print(df.format(0)+", ");
109             }
110             System.out.print(df.format(1));
111             for (int i=j+1;i<n;i++) {
112                 System.out.print(", "+df.format(0));
113             }
114         } else {
115             double len = 1;
116             for (int i=0;i<r;i++) {
117                 len += mat[i][j]*mat[i][j];
118             }
119             len = Math.sqrt(len);
120             System.out.print(df.format(mat[0][j]/len));
121             for (int i=1;i<r;i++) {
122                 System.out.print(", ");
123                 System.out.print(df.format(mat[i][j]/len));
124             }
125             for (int i=r;i<j;i++) {
126                 System.out.print(", "+df.format(0));
127             }
128             System.out.print(", "+df.format(-1/len));
129             for (int i=j+1;i<n;i++) {
130                 System.out.print(", "+df.format(0));
131             }
132         }
133         System.out.println(")‘");
134     }
135     private boolean zero(double x) {
136         return x<0.000001 && x>-0.000001;
137     }
138 }
139 
140 public class LinearEquation {
141     
142     public static void main(String[] args) {
143         Scanner in = new Scanner(System.in);
144         System.out.print("Number of Equations: ");
145         int m = in.nextInt();
146         System.out.print("Number of Variables: ");
147         int n = in.nextInt();
148         Linear eq = new Linear(m,n);
149         System.out.println("Augmented Matrix:");
150         eq.getData(in);
151         in.close();
152         eq.solve();
153     }
154 }

 

2. Orthogonal Diagonalization

  1 import java.util.*;
  2 import java.text.*;
  3 
  4 class Vector {
  5     private int dim;
  6     private double[] data;
  7     
  8     public Vector(int dim) {
  9         this.dim = dim;
 10         data = new double[dim];
 11     }
 12     public Vector(Vector other) {
 13         dim = other.dim;
 14         data = new double[dim];
 15         for (int i=0;i<dim;i++) {
 16             data[i] = other.data[i];
 17         }
 18     }
 19     public int getDim() {
 20         return dim;
 21     }
 22     public double get(int i) {
 23         if (i<0||i>=dim) {
 24             throw new RuntimeException("Index out of Bound");
 25         }
 26         return data[i];
 27     }
 28     public void set(int i,double k) {
 29         if (i<0||i>=dim) {
 30             throw new RuntimeException("Index out of Bound");
 31         }
 32         data[i] = k;
 33     }
 34     public Vector add(Vector other) {
 35         if (other.dim!=dim) {
 36             throw new RuntimeException("Invalid Vector Addition");
 37         }
 38         Vector val = new Vector(dim);
 39         for (int i=0;i<dim;i++) {
 40             val.data[i] = data[i]+other.data[i];
 41         }
 42         return val;
 43     }
 44     public Vector mult(double k) {
 45         Vector val = new Vector(dim);
 46         for (int i=0;i<dim;i++) {
 47             val.data[i] = data[i]*k;
 48         }
 49         return val;
 50     }
 51     public double mult(Vector other) {
 52         if (other.dim!=dim) {
 53             throw new RuntimeException("Invalid Vector Multiplication");
 54         }
 55         double val = 0;
 56         for (int i=0;i<dim;i++) {
 57             val += data[i]*other.data[i];
 58         }
 59         return val;
 60     }
 61     public Vector unify() {
 62         double len = 0;
 63         for (int i=0;i<dim;i++) {
 64             len += data[i]*data[i];
 65         }
 66         len = Math.sqrt(len);
 67         Vector val = new Vector(this);
 68         if (zero(len)) {
 69             return val;
 70         } else if (data[0]<0) {
 71             len = -len;
 72         }
 73         for (int i=0;i<dim;i++) {
 74             val.data[i] /= len;
 75         }
 76         return val;
 77     }
 78     private boolean zero(double k) {
 79         return k>-0.000001&&k<0.000001;
 80     }
 81 }
 82 
 83 class Matrix {
 84     private int size;
 85     private double[][] data;
 86     
 87     public Matrix(int size) {
 88         this.size = size;
 89         data = new double[size][size];
 90     }
 91     public Matrix(Matrix other) {
 92         size = other.size;
 93         data = new double[size][size];
 94         for (int i=0;i<size;i++) {
 95             for (int j=0;j<size;j++) {
 96                 data[i][j] = other.data[i][j];
 97             }
 98         }
 99     }
100     public int getSize() {
101         return size;
102     }
103     public double get(int i,int j) {
104         if (i<0||i>=size||j<0||j>=size) {
105             throw new RuntimeException("Index out of Bound");
106         }
107         return data[i][j];
108     }
109     public void set(int i,int j,double k) {
110         if (i<0||i>=size||j<0||j>=size) {
111             throw new RuntimeException("Index out of Bound");
112         }
113         data[i][j] = k;
114     }
115     public double trace() {
116         double val = 0;
117         for (int i=0;i<size;i++) {
118             val += data[i][i];
119         }
120         return val;
121     }
122     public double det() {
123         double val = 0;
124         if (size==1) {
125             val += data[0][0];
126         } else if (size==2) {
127             val += data[0][0]*data[1][1];
128             val -= data[0][1]*data[1][0];
129         } else if (size==3){
130             val += data[0][0]*data[1][1]*data[2][2];
131             val += data[0][1]*data[1][2]*data[2][0];
132             val += data[1][0]*data[2][1]*data[0][2];
133             val -= data[0][2]*data[1][1]*data[2][0];
134             val -= data[1][2]*data[2][1]*data[0][0];
135             val -= data[2][2]*data[0][1]*data[1][0];
136         } else {
137             throw new RuntimeException("High-Dim not Supported");
138         }
139         return val;
140     }
141     public void getEig(double[] lambda,Vector[] q) {
142         // Get Eigenvalues and Eigenvetors:
143         if (size==1) {
144             lambda[0] = data[0][0];
145         } else if (size==2) {
146             equSolve(lambda,-trace(),det());
147         } else if (size==3) {
148             double c = get(0,0)*get(1,1)+get(1,1)*get(2,2)+get(2,2)*get(0,0)
149                     - get(0,1)*get(1,0) - get(1,2)*get(2,1) - get(2,0)*get(0,2);
150             equSolve(lambda,-trace(),c,-det());
151             if (zero(lambda[0]-lambda[2])) {
152                 double tmp = lambda[2];
153                 lambda[2] = lambda[1];
154                 lambda[1] = tmp;
155             }
156         }
157         for (int i=0;i<size;i+=getEigvect(lambda,q,i));
158         schmidt(q);
159     }
160     private int getEigvect(double[] lambda,Vector[] q,int pos) {
161         int val = 1;
162         for (;pos+val<size;val++) {
163             if (!zero(lambda[pos+val]-lambda[pos])) {
164                 break;
165             }
166         }
167         double[][] mat = new double[size][size];
168         for (int i=0;i<size;i++) {
169             for (int j=0;j<size;j++) {
170                 mat[i][j] = -data[i][j];
171             }
172             mat[i][i] += lambda[pos];
173         }
174         for (int j=0;j<size;j++) {
175             utilize(mat,j);
176             for (int i=j+1;i<size;i++) {
177                 subtract(mat,i,j);
178             }
179         }
180         int r = size-val;
181         for (int i=r-2;i>=0;i--) {
182             for (int j=i+1;j<r;j++) {
183                 subtract(mat,i,j);
184             }
185         }
186         for (int j=r;j<size;j++) {
187             for (int i=0;i<r;i++) {
188                 q[pos+j-r].set(i,mat[i][j]);
189             }
190             q[pos+j-r].set(j,-1);
191         }
192         return val;
193     }
194     private void equSolve(double[] rt,double b,double c) {
195         // Try to solve a Quadratic Equation: x^2+b*x+c = 0
196         double delt = b*b-4*c;
197         if (delt<0) {
198             throw new RuntimeException("No Real Solution");
199         }
200         delt = Math.sqrt(delt);
201         rt[0] = (-b-delt)/2;
202         rt[1] = (-b+delt)/2;
203     }
204     private void equSolve(double[] rt,double b,double c,double d) {
205         // Try to solve a Cubic Equation:  x^3+b*x^2+c*x+d = 0
206         double alpha = -b*b*b/27-d/2+b*c/6;
207         double beta = c/3-b*b/9;
208         if (alpha*alpha+beta*beta*beta>0) {
209             throw new RuntimeException("No Real Solution");
210         } else if (zero(beta)){
211             rt[0] = rt[1] = rt[2] = -b/3; 
212             return;
213         }
214         double tmp = Math.acos(alpha/Math.sqrt(-beta*beta*beta));
215         rt[0] = -b/3+2*Math.sqrt(-beta)*Math.cos((tmp-2*Math.PI)/3);
216         rt[1] = -b/3+2*Math.sqrt(-beta)*Math.cos(tmp/3);
217         rt[2] = -b/3+2*Math.sqrt(-beta)*Math.cos((tmp+2*Math.PI)/3);
218     }
219     private void schmidt(Vector[] vect) {
220         // Gram–Schmidt Orthogonalization
221         int num = vect.length;
222         Vector[] tmp = new Vector[num];
223         for (int i=0;i<num;i++) {
224             tmp[i] = new Vector(vect[i]);
225         }
226         for (int i=0;i<num;i++) {
227             for (int j=0;j<i;j++) {
228                 double k = tmp[i].mult(vect[j]);
229                 k /= vect[j].mult(vect[j]);
230                 vect[i] = vect[i].add(vect[j].mult(-k));
231             }
232             vect[i] = vect[i].unify();
233         }
234     }
235     private void utilize(double[][] mat,int j) {
236         int t = size-1;
237         for (int i=j;i<=t;i++) {
238             if (zero(mat[i][j])) {
239                 for (int k=j;k<size;k++) {
240                     double tmp = mat[i][k];
241                     mat[i][k] = mat[t][k];
242                     mat[t][k] = tmp;
243                 }    i --;    t --;
244             } else {
245                 for (int k=size-1;k>=j;k--) {
246                     mat[i][k]/=mat[i][j];
247                 }
248             }
249         }
250     }
251     private void subtract(double[][] mat,int i,int j) {
252         if (zero(mat[j][j])) {
253             return;
254         }
255         for (int k=size-1;k>=j;k--) {
256             mat[i][k] -= mat[i][j]*mat[j][k];
257         }
258     }
259     private boolean zero(double k) {
260         return k>-0.000001&&k<0.000001;
261     }
262 }
263 
264 public class Main {
265     public static Scanner in;
266     public static DecimalFormat df;
267     
268     public static void getInput(Matrix coef,Vector rear) {
269         int dim = coef.getSize();
270         System.out.println("Get the Equation:");
271         System.out.print("0 = \t");
272         // Second-Degree Terms:
273         System.out.print("x^2 * ");
274         coef.set(0,0,in.nextDouble());
275         System.out.print("  +\t");
276         System.out.print("y^2 * ");
277         coef.set(1,1,in.nextDouble());
278         System.out.print("  +\t");
279         if (dim==3) {
280             System.out.print("z^2 * ");
281             coef.set(2,2,in.nextDouble());
282             System.out.print("  +\t");
283         }
284         System.out.print("x*y * ");
285         coef.set(0,1,in.nextDouble()/2);
286         coef.set(1,0,coef.get(0,1));
287         System.out.print("  +\t");
288         if (dim==3) {
289             System.out.print("y*z * ");
290             coef.set(1,2,in.nextDouble()/2);
291             coef.set(2,1,coef.get(1,2));
292             System.out.print("  +\t");
293             System.out.print("z*x * ");
294             coef.set(2,0,in.nextDouble()/2);
295             coef.set(0,2,coef.get(2,0));
296             System.out.print("  +\t");
297         }
298         // First-Degree Terms:
299         System.out.print("x * ");
300         rear.set(0,in.nextDouble());
301         System.out.print("  +\t");
302         System.out.print("y * ");
303         rear.set(1,in.nextDouble());
304         System.out.print("  +\t");
305         if (dim==3) {
306             System.out.print("z * ");
307             rear.set(2,in.nextDouble());
308             System.out.print("  +\t");
309         }
310         rear.set(dim,in.nextDouble());
311     }
312     public static void solve(Matrix coef,Vector rear) {
313         int dim = coef.getSize();
314         double[] lambda = new double[dim];
315         Vector[] q = new Vector[dim];
316         for (int i=0;i<dim;i++) {
317             q[i] = new Vector(dim);
318         }
319         coef.getEig(lambda,q);
320         printRes(lambda,q,rear);
321     }
322     public static void printRes(double[] lambda,Vector[] q,Vector rear) {
323         int dim = lambda.length;
324         double[] bias = new double[dim];
325         double cst = rear.get(dim);
326         for (int i=0;i<dim;i++) {
327             for (int j=0;j<dim;j++) {
328                 bias[i]+=q[i].get(j)*rear.get(j);
329             }
330             bias[i] /= (2*lambda[i]);
331             cst -= lambda[i]*bias[i]*bias[i];
332         }
333         System.out.print("0 = \t");
334         System.out.print(df.format(lambda[0]));
335         System.out.print(" * (");
336         translate(0,q,bias[0]);
337         System.out.println(")^2");
338         System.out.print("  +\t");
339         System.out.print(df.format(lambda[1]));
340         System.out.print(" * (");
341         translate(1,q,bias[1]);
342         System.out.println(")^2");
343         if (dim==3) {
344             System.out.print("  +\t");
345             System.out.print(df.format(lambda[2]));
346             System.out.print(" * (");
347             translate(2,q,bias[2]);
348             System.out.println(")^2");
349         }
350         System.out.println("  +\t"+df.format(cst));
351     }
352     private static void translate(int pos,Vector[] q,double c) {
353         System.out.print(df.format(q[pos].get(0))+"*x");
354         System.out.print("+"+df.format(q[pos].get(1))+"*y");
355         if (q.length==3) {
356             System.out.print("+"+df.format(q[pos].get(2))+"*z");
357         }
358         if (!zero(c)) {
359             System.out.print("+"+df.format(c));
360         }
361     }
362     private static boolean zero(double k) {
363         return k>-0.000001&&k<0.000001;
364     }
365     public static void main(String[] args) {
366         in = new Scanner(System.in);
367         df = new DecimalFormat("#.0000");
368         System.out.println("Curve or Surface?");
369         System.out.print("\tdim = ");
370         int dim = in.nextInt();
371         if (dim>=2 && dim<=3) {
372             Matrix coef = new Matrix(dim);
373             Vector rear = new Vector(dim+1);
374             getInput(coef,rear);
375             solve(coef,rear);
376         } else {
377             System.out.println("Sorry, I can‘t do that.");
378         }
379         in.close();
380     }
381 }

 

 

======================== 不相亏欠,何须怀缅 =========================

 

  人只有将寂寞坐断

  才可以重拾喧闹

  把悲伤过尽

  才可以重见欢颜

  把苦涩尝遍

  就会自然回甘

        —— 林徽因

Commemoration of Linear Algebra

标签:

原文地址:http://www.cnblogs.com/DevinZ/p/4435859.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!