标签:
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