From 13c23aafa89bd19470f0d109ea8447fb48322ff6 Mon Sep 17 00:00:00 2001 From: Daniel Weschke Date: Sat, 15 Mar 2014 10:47:42 +0100 Subject: [PATCH] Add some polynomial vector and number classs --- src/awt/Draw.java | 593 ++++++++ src/math/equation/Polynomial.java | 1195 +++++++++++++++++ src/math/equation/PolynomialComplex.java | 114 ++ src/math/equation/RationalPolynomial.java | 128 ++ .../equation/RationalPolynomialComplex.java | 112 ++ src/math/matrix/VectorComplex.java | 513 +++++++ src/math/number/Number.java | 12 + src/math/number/NumberComplex.java | 377 ++++++ src/physics/Nyquist.java | 154 +++ 9 files changed, 3198 insertions(+) create mode 100644 src/awt/Draw.java create mode 100644 src/math/equation/Polynomial.java create mode 100644 src/math/equation/PolynomialComplex.java create mode 100644 src/math/equation/RationalPolynomial.java create mode 100644 src/math/equation/RationalPolynomialComplex.java create mode 100644 src/math/matrix/VectorComplex.java create mode 100644 src/math/number/Number.java create mode 100644 src/math/number/NumberComplex.java create mode 100644 src/physics/Nyquist.java diff --git a/src/awt/Draw.java b/src/awt/Draw.java new file mode 100644 index 0000000..cbdd353 --- /dev/null +++ b/src/awt/Draw.java @@ -0,0 +1,593 @@ +package awt; + +import javax.swing.JLabel; +import javax.swing.SwingConstants; +import javax.swing.SwingUtilities; + +import math.Maths; +import math.matrix.Matrix; +import math.matrix.Vector; +import physics.Nyquist; +import stdlib.Color; +import stdlib.StdDraw; + +public class Draw extends StdDraw{ + static Thread plot; + static JLabel infoLabel; + + public static double scaleLeft = -10; + public static double scaleRight = 10; + public static double scaleBottom = -10; + public static double scaleTop = 10; + + // TODO: combine X,Y,XM,YM ??? + // see in animate + static double[] X; + static double[] Y; + static Matrix XM; + static Matrix YM; + + static Matrix M; + static Matrix dM; + static Vector m; // for FEM + static double fac; // for FEM + + boolean animationSet = false; + + private static double process; + + /** + * plot + */ + public Draw(){ + setupGUI(); + } + + /** + * plot + */ + public Draw(String title){ + setupGUI(title); + } + + /** + * xy-plot + * @param X + * @param Y + */ + public Draw(int[] X, int[] Y){ + setupGUI(); + allocation(X, Y); + setScale(Maths.min(X)*1.1, Maths.max(X)*1.1, Maths.min(Y)*1.1, Maths.max(Y)*1.1); + } + + /** + * xy-plot + * @param X + * @param Y + */ + public Draw(int[] X, int[] Y, String title){ + setupGUI(title); + allocation(X, Y); + setScale(Maths.min(X)*1.1, Maths.max(X)*1.1, Maths.min(Y)*1.1, Maths.max(Y)*1.1); + } + + /** + * xy-plot + * @param X + * @param Y + */ + public Draw(double[] X, double[] Y){ + setupGUI(); + allocation(X, Y); + setScale(Maths.min(X)*1.1, Maths.max(X)*1.1, Maths.min(Y)*1.1, Maths.max(Y)*1.1); + } + + /** + * xy-plot + * @param X + * @param Y + */ + public Draw(double[] X, double[] Y, String title){ + setupGUI(title); + allocation(X, Y); + setScale(Maths.min(X)*1.1, Maths.max(X)*1.1, Maths.min(Y)*1.1, Maths.max(Y)*1.1); + } + + /** + * xy-plot + * @param X + * @param Y + */ + public Draw(Vector X, Vector Y){ + setupGUI(); + allocation(X, Y); + setScale(X.min()*1.1, X.max()*1.1, Y.min()*1.1, Y.max()*1.1); + } + + /** + * xy-plot + * @param X + * @param Y + */ + public Draw(Vector X, Vector Y, String title){ + setupGUI(title); + allocation(X, Y); + setScale(X.min()*1.1, X.max()*1.1, Y.min()*1.1, Y.max()*1.1); + } + + /** + * xy-plot + * @param X + * @param Y + */ + public Draw(Matrix X, Matrix Y){ + setupGUI(); + allocation(X, Y); + setScale(X.min()*1.1, X.max()*1.1, Y.min()*1.1, Y.max()*1.1); + } + + /** + * xy-plot + * @param X + * @param Y + */ + public Draw(Matrix X, Matrix Y, String title){ + setupGUI(title); + allocation(X, Y); + setScale(X.min()*1.1, X.max()*1.1, Y.min()*1.1, Y.max()*1.1); + } + + +// GUI + + private void setupGUI(){ + frame.setLocation(100, 50); + } + + private void setupGUI(String title){ + frame.setLocation(100, 50); + frame.setTitle(title); + } + + public void setBGC(Color bgColor){ + clear(bgColor); + } + + public void setFGC(Color Color){ + setGraphColor(Color); + } + + +// Graph + + private void allocation(int[] X, int[] Y){ + int imax = X.length; + Draw.XM = new Matrix(imax,1); + Draw.YM = new Matrix(imax,1); + int i; + for(i=0; i class: DrawFEM? + @SuppressWarnings("deprecation") + public void animateFEM(double[][] A, double[][] dA, double[] a) { + M = new Matrix(A); + dM = new Matrix(dA); + m = new Vector(a); + // Display size + double xMin = M.getN(2).min(); + double xMax = M.getN(2).max(); + double yMin = M.getN(3).min(); + double yMax = M.getN(3).max(); + if(M.getN(0).min() < xMin) xMin = M.getN(0).min(); + if(M.getN(0).max() > xMax) xMax = M.getN(0).max(); + if(M.getN(1).min() < yMin) yMin = M.getN(1).min(); + if(M.getN(1).max() > yMax) yMax = M.getN(1).max(); +// System.out.println((xMax-xMin)+" "+(yMax-yMin)); +// System.out.println(((xMax-xMin) - (yMax-yMin))/2); +// System.out.println(yMin+" "+yMax); + double xMM = xMax - xMin; + double yMM = yMax - yMin; + if((xMM) > (yMM)){ + double dif = (xMM - yMM)/2; + yMin -= dif*1.05; // kleine Verschiebung nach oben + yMax += dif*0.95; + } else { + double dif = (yMM - xMM)/2; + xMin -= dif*1.05; + xMax += dif*0.95; + } +// System.out.println(yMin+" "+yMax); +// System.out.println((xMax-xMin)+" "+(yMax-yMin)); + setScale(xMin, xMax, yMin, yMax); + // Displacement factor + double xdMM = dM.getN(2).max() - dM.getN(2).min(); + double ydMM = dM.getN(3).max() - dM.getN(3).min(); + double dMM = Math.sqrt(xdMM*xdMM+ydMM*ydMM); + double MM = Math.sqrt(xMM*xMM+yMM*yMM); + fac = MM/dMM/100*2; + if(fac<1)fac=1; +// System.out.println(dMM); +// System.out.println(MM); + System.out.printf("Scalierung %.4g%n",fac); + if(!animationSet){ + animateContents(); + SwingUtilities.invokeLater(new Runnable() { + @Override + public void run() { + setPenColor(Color.LIGHT_GRAY); + new AnimateFEM(); + } + }); + animationSet = true; + } else contents(); + setPenColor(Color.LIGHT_GRAY); + if(plot!=null) plot.stop(); + plot = new Thread(new AnimateFEM()); + plot.start(); + } + + private void contents(){ + frame.setLocation(100, 50); + + setXscale(scaleLeft, scaleRight); + setYscale(scaleBottom, scaleTop); + } + + public void cross(){ + setPenColor(Color.GRAY); + line(scaleLeft, 0, scaleRight, 0); + line(0, scaleBottom, 0, scaleTop); + setPenColor(getGraphColor()); + } + + private void animateContents(){ + contents(); + infoLabel = new JLabel("0 %"); + infoLabel.setHorizontalAlignment(SwingConstants.RIGHT); + infoLabel.setBounds(DEFAULT_SIZE+left-107, top+5, 100, 20); + frame.add(infoLabel); + } + + + /** + * Initialize Nyquist data + */ + private void nyquist(){ + frame.setTitle("Nyquist"); + Nyquist nyq = new Nyquist(); + Draw.X = nyq.getNyquist().getN(0).getArray(); + Draw.Y = nyq.getNyquist().getN(1).getArray(); + + scaleLeft = -1.5; + scaleRight = 2.5; + scaleBottom = -2.5; + scaleTop = 1.5; + } + + + /** + * Plot Nyquist + */ + public void plotNyquist() { + nyquist(); + plot(); + } + + + /** + * Animate Nyquist + */ + public void animateNyquist() { + nyquist(); + animate(); + cross(); + //TODO:mark? + markX(-1); + } + + /** + * + */ + public class Animate implements Runnable{ + @Override + public void run(){ + if(XM != null && YM != null){ + // TODO: different Color + Color[] color = new Color[]{ + Color.LIGHT_BLUE, + Color.PURPLE, + Color.YELLOW, + Color.GREEN, + Color.RED}; + int i, j; + int imax = XM.getM(); + int jmax = XM.getN(); + int max = imax > jmax ? imax : jmax; + int sec = (int) 1000./max*10; + sec = sec > 10 ? 10 : sec; +// System.out.println(max); +// System.out.println(sec); + for(j=0; j=color.length) setPenColor(color[0]); + else setPenColor(color[i]); + point(XM.get(i,j), YM.get(i,j)); + show(sec); + setProcess((double)(i+j+1)/(imax+jmax)*100); + } + } + } else if(X!=null){ + int i; + int imax = X.length; + int sec = (int) 1000./imax*10; + sec = sec > 10 ? 10 : sec; +// System.out.println(imax); +// System.out.println(sec); + for(i=0; i class: DrawFEM? + public class AnimateFEM implements Runnable{ + @Override + public void run(){ + if(M!=null){ + int i, j; + int imax = M.getM(); + int jmax = 200; + int mMax = m.abs().indexOfMax(); + double max = m.max(); + double min = m.min(); + double f; + for(j=0; j= 0.92*max) setPenColor(Color.ORANGE); + else if(m.get(i) >= 0.81*max) setPenColor(Color.YELLOW); + else if(m.get(i) >= 0.67*max) setPenColor(Color.GREEN); + else if(m.get(i) >= 0.50*max) setPenColor(Color.GREEN); + else if(m.get(i) >= 0.30*max) setPenColor(Color.GREEN); + else if(m.get(i) > 0.30*min) setPenColor(Color.GREEN); + else if(m.get(i) > 0.50*min) setPenColor(Color.GREEN); + else if(m.get(i) > 0.67*min) setPenColor(Color.GREEN); + else if(m.get(i) > 0.81*min) setPenColor(Color.GREEN); + else if(m.get(i) > 0.92*min) setPenColor(Color.CYAN); + else if(m.get(i) > 1.00*min) setPenColor(Color.BLUE); + else setPenColor(Color.PURPLE); + line(M.get(i,0)+dM.get(i,0)*f, M.get(i,1)+dM.get(i,1)*f, + M.get(i,2)+dM.get(i,2)*f, M.get(i,3)+dM.get(i,3)*f); + } + show(10); + setProcess((double)(i+j+1)/(imax+jmax)*100); + } + } + } + } + + @SuppressWarnings("deprecation") + public void stop(){ + plot.stop(); +// plot.interrupt(); +// frame.repaint(); // bringt nichts +// frame.removeAll(); + } + + public void markX(double x){ + double size = (scaleTop-scaleBottom)*.01; + line(x, -size, x, size); + if(x!=0) + if(scaleLeft==0) text(x, 3*size, String.format("%.1f", x)); + else text(x, -4*size, String.format("%.1f", x)); + } + + public void markY(double y){ + double size = (scaleRight-scaleLeft)*.01; + line(-size, y, size, y); + if(y!=0) + if(scaleBottom==0) textLeft(2*size, y, String.format("%.1f", y)); + else textRight(-2.5*size, y, String.format("%.1f", y)); + } + + public void mark(double x, double y){ + double sizeX = (scaleTop-scaleBottom)*.01; + double sizeY = (scaleRight-scaleLeft)*.01; + filledCircle(x, y, sizeXi , ordered by + * ann + an-1n-1 + ... + a2x2 + a1x + a0 + */ + private double[] coef; + + /** + * integration constants ci ordered by + * cnn + cn-1n-1 + ... + c2x2 + c1x + c0 + */ + private double[] constI; + + /** + * Construct a zero polynomial, 0. + */ + public Polynomial(){ + this(0,0); + } + + /** + * Construct a zero polynomial, + * 0x0 + 0x1 + ... + 0xn. + * @param n number of zero coefficient + */ + public Polynomial(int n){ + coef = new double[n]; + } + + /** + * Construct a constant polynomial, + * ax0. + * @param a constant + */ + public Polynomial(double a){ + coef = new double[]{a}; + } + + /** + * Construct a polynomial, + * axb. + * @param a coefficient + * @param b exponent + */ + public Polynomial(double a, int b) { + coef = new double[b+1]; + coef[b] = a; + } + + /** + * Construct a linear polynomial, + * a1x + a0. + * @param a1 coefficient + * @param a0 coefficient + */ + public Polynomial(double a1, double a0) { + coef = new double[]{a0, a1}; + } + + /** + * Construct a quadratic polynomial, + * a2x2 + a1x + a0. + * @param a2 coefficient + * @param a1 coefficient + * @param a0 coefficient + */ + public Polynomial(double a2, double a1, double a0) { + coef = new double[]{a0, a1, a2}; + } + + /** + * Construct a cubic polynomial, + * a3x3 + a2x2 + a1x + a0. + * @param a3 coefficient + * @param a2 coefficient + * @param a1 coefficient + * @param a0 coefficient + */ + public Polynomial(double a3, double a2, double a1, double a0) { + coef = new double[]{a0, a1, a2, a3}; + } + + /** + * Construct a polynomial with an array of coefficients. + * Ordered by: ann + an-1n-1 + ... + a2x2 + a1x + a0 + * @param a coefficients + */ + public Polynomial(double[] a){ + this(a, true); + } + + /** + * Construct a polynomial with an array of coefficients. + * Ordered by: ann + an-1n-1 + ... + a2x2 + a1x + a0 + * @param a coefficients + * @param reverse order + */ + public Polynomial(double[] a, boolean reverse){ + int size = a.length; + coef = new double[size]; + int i; + if(reverse) + for(i=0; inn + an-1n-1 + ... + a2x2 + a1x + a0 + * @param a coefficients + */ + public Polynomial(Object... a){ + this(true, a); + } + + /** + * Construct a polynomial with an array of coefficients. + * Ordered by: ann + an-1n-1 + ... + a2x2 + a1x + a0 + * @param a coefficients + */ + public Polynomial(boolean reverse, Object[] a){ + int size = a.length; + coef = new double[size]; + int i; + if(reverse) + for(i=0; inn + an-1n-1 + ... + a2x2 + a1x + a0 + * @param a coefficients + */ + public Polynomial(String a){ + this(0,0); + int i; + + String[] p = a.trim().split("\\s+"); + int size = p.length; + + if(size > 0){ + coef = new double[size]; + for(i=0; i + * f(px) = a = py. + * @param p + */ + public Polynomial(Point2D p){ + coef = new double[]{p.y()}; + } + + /** + * Linear function. + * f(x) = a0 + a1x .
+ * f(p1x) = a0 + a1p1x = p1y
+ * f(p2x) = a0 + a1p2x = p2y + * @param p1 + * @param p2 + */ + public Polynomial(Point2D p1, Point2D p2){ + Matrix A = new Matrix(new double[][]{{1, p1.x()}, {1, p2.x()}}); + Vector b = new Vector(new double[]{p1.y(), p2.y()}); + coef = A.solve(b).getArray(); + } + + public Polynomial(Point2D p, double alpha){ + this(p, alpha, false); + } + + public Polynomial(Point2D p, double alpha, boolean degree){ + double a = alpha; + if(degree) a = Maths.cosd(alpha); + coef = new double[]{p.y()-a*p.x(), a}; + } + + /** + * Copy constructor + * @param a polynomial + */ + public Polynomial(Polynomial a){ + coef = a.coef; + constI = a.constI; + } + + /** + * Binomial (1+x)n + *

+ * n = 0 : 1x^0
+ * n = 3 : 1x^3 + 3x^2 + 3x^1 + 1x^0
+ * n = 6 : 1x^6 + 6x^5 + 15x^4 + 20x^3 + 15x^2 + 6x^1 + 1x^0
+ * @param n + * @return (1+x)n + */ + public static Polynomial binomial(int n){ + Polynomial one = new Polynomial(1, 0); // 1 + Polynomial x = new Polynomial(1, 1); // x + Polynomial binomial = one; // 1 + int i; + for(i=0; i + * T0 = 1
+ * T1 = x
+ * Tn = 2x * Tn-1 - Tn-2
+ *

+ * n = 1: 1x0
+ * n = 2: 1x1
+ * n = 9: 256x9 - 576x7 + 432x5 - 120x3 + 9x1
+ * @param n th Chebyshev polynomial + * @return Tn = 2x * Tn-1 - Tn-2 + * @see Polynomials#chebyshevT(int) + */ + public static Polynomial chebyshevT(int n){ + if(n == 0) return new Polynomial(1, 0); // 1 + else if(n == 1) return new Polynomial(1, 1); // x + else return new Polynomial(2, 1).times(chebyshevT(n-1)).minus(chebyshevT(n-2)); + } + + /** + * n-th Chebyshev polynomial of the second kind. + *

+ * U0 = 1
+ * U1 = 2x
+ * Un = 2x * Un-1 - Un-2
+ *

+ * n = 1: 1x0
+ * n = 2: 2x1
+ * n = 9: 512x9 - 1024x7 + 672x5 - 160x3 + 10x1
+ * @param n th Chebyshev polynomial + * @return Un = 2x * Un-1 - Un-2 + * @see Polynomials#chebyshevU(int) + */ + public static Polynomial chebyshevU(int n){ + if(n == 0) return new Polynomial(1, 0); // 1 + else if(n == 1) return new Polynomial(2, 1); // 2x + else return new Polynomial(2, 1).times(chebyshevU(n-1)).minus(chebyshevU(n-2)); + } + + /** + * n-th Hermite polynomial. + *

+ * H0 = 1
+ * H1 = 2x
+ * Hn = 2x * Hn-1 - 2(n-1) * Hn-2
+ *

+ * n = 1: 1x0
+ * n = 2: 2x1
+ * n = 9: 512x9 - 9216x7 + 48384x5 - 80640x3 + 30240x1
+ * @param n th Hermite polynomial + * @return Hn = 2x * Hn-1 - 2(n-1) * Hn-2 + * @see Polynomials#hermite(int) + */ + public static Polynomial hermite(int n){ + if(n == 0) return new Polynomial(1, 0); // 1 + else if(n == 1) return new Polynomial(2, 1); // 2x + else return new Polynomial(2, 1).times(hermite(n-1)).minus(hermite(n-2).times(2*(n-1))); + } + + /** + * First derivative of n-th Hermite polynomial. + *

+ * H'0 = 0
+ * H'n = 2n * Hn-1
+ *

+ * n = 1: 2x0
+ * n = 2: 8x1
+ * n = 9: 4608x8 - 64512x6 + 241920x4 - 241920x2 + 30240x0
+ * @param n th Hermite polynomial + * @return H'n = 2n * Hn-1 + * @see #hermiteD(int, int) + */ + public static Polynomial hermiteD(int n){ + if(n == 0) return new Polynomial(0, 0); // 0 + else return hermite(n-1).times(2*n); + } + + /** + * m-th derivative of n-th Hermite polynomial. + *

+ * n = 9, m = 2: 36864x7 - 387072x5 + 967680x3 - 483840x1
+ * @param n th Hermite polynomial + * @param m derivative + * @return H(m)n = 2n * H(m-1)n-1 + * @see #hermiteD(int) + */ + public static Polynomial hermiteD(int n, int m){ + if(m == 0) return hermite(n); + if(n < m) return new Polynomial(0, 0); // 0 + if(m == 1) return hermiteD(n); + else return hermiteD(n-1,m-1).times(2*n); + } + + /** + * n-th Hermite polynomial (probabilist). + *

+ * He0 = 1
+ * He1 = x
+ * Hen = x * Hen-1 - (n-1) * Hen-2
+ *

+ * n = 1: 1x0
+ * n = 2: x1
+ * n = 9: x9 - 36x7 + 378x5 - 1260x3 + 945x1
+ * @param n th Hermite polynomial + * @return Hen = x * Hen-1 - (n-1) * Hen-2 + * @see Polynomials#hermiteE(int) + */ + public static Polynomial hermiteE(int n){ + if(n == 0) return new Polynomial(1, 0); // 1 + else if(n == 1) return new Polynomial(1, 1); // x + else return new Polynomial(1, 1).times(hermiteE(n-1)).minus(hermiteE(n-2).times(n-1)); + } + + /** + * First derivative of n-th Hermite polynomial (probabilist). + *

+ * He'0 = 0
+ * He'n = n * Hen-1
+ *

+ * n = 1: 1x0
+ * n = 2: 2x1
+ * n = 9: 9x8 - 252x6 + 1890x4 - 3780x2 + 945x0
+ * @param n th Hermite polynomial + * @return He'n = n * Hen-1 + * @see #hermiteED(int, int) + */ + public static Polynomial hermiteED(int n){ + if(n == 0) return new Polynomial(0, 0); // 0 + else return hermiteE(n-1).times(n); + } + + /** + * m-th derivative of n-th Hermite polynomial (probabilist). + *

+ * n = 9, m = 2: 72x7 - 1512x5 + 7560x3 - 7560x1
+ * @param n th Hermite polynomial + * @param m derivative + * @return He(m)n = n * He(m-1)n-1 + * @see #hermiteED(int) + */ + public static Polynomial hermiteED(int n, int m){ + if(m == 0) return hermiteE(n); + if(n < m) return new Polynomial(0, 0); // 0 + if(m == 1) return hermiteED(n); + else return hermiteED(n-1,m-1).times(n); + } + + /** + * + * @return the degree of this polynomial (0 for the zero polynomial) + */ + public int degree(){ + int d = 0, i; + if(coef != null) + for(i=0; i=0; i--) { + Polynomial term = new Polynomial(a.coef[i], 0); + c = term.plus(b.times(c)); + } + return c; + } + + /** + * use Horner's method to compute and return the polynomial evaluated at x + * @param x + * @return the polynomial evaluated at x + */ + public double evaluate(double x){ + int i; + double p = 0; + for(i=degree(); i>=0; i--) + p = coef[i] + (x * p); + return p; + } + + /** + * Compute the polynomial evaluation at xi using Horner's method. + * @param x + * @return the polynomial evaluation at xi + */ + public double[] evaluate(double[] x){ + int i; + double[] result = new double[x.length]; + for(i=0;ii using Horner's method. + * @param x + * @return the polynomial evaluation at xi + * @throws IllegalDimensionException x must no be null + */ + public Vector evaluate(Vector x) throws IllegalDimensionException{ + if(x == null) throw new IllegalDimensionException(); + int i; + Vector result = new Vector(x.n()); + for(i=0;ii using Horner's method. + * @param start + * @param increment + * @param end + * @return the polynomial evaluation at xi + * @throws IllegalDimensionException x must no be null + */ + public Matrix2D evaluateAsMatrix2D(double start, double increment, double end) throws IllegalDimensionException{ + Vector x = Vector.fill(start, increment, end); + return evaluateAsMatrix2D(x); + } + + /** + * Compute the polynomial evaluation at xi using Horner's method. + * @param x + * @return the polynomial evaluation at xi + * @throws IllegalDimensionException x must no be null + */ + public Matrix2D evaluateAsMatrix2D(Vector x) throws IllegalDimensionException{ + if(x == null) throw new IllegalDimensionException(); + int i; + Matrix2D result = new Matrix2D(x.n()); +// System.out.println(x.length()); +// System.out.println(result.getM() + " * " + result.getN()); + for(i=0;i0) deriv = derive().derive(n-1); + return deriv; + } + + /** + * derive this polynomial + * @param x + * @return f'(x) + */ + public double derive(double x){ + return derive().evaluate(x); + } + + /** + * derive this polynomial n times + * @param n + * @param x + * @return f^n(x) + */ + public double derive(int n, double x){ + return derive(n).evaluate(x); + } + + /** + * integrate this polynomial and return it + * @return F(x) + */ + public Polynomial integrate(){ + Polynomial integ; + if(degree() == 0 && getCoef(0) == 0){ + integ = new Polynomial(coef.length+1); + } else { + integ = new Polynomial(0, degree() + 1); + for(int i=0; i<=degree(); i++){ + integ.coef[i+1] = 1./(i+1) * coef[i]; + // System.out.println("1/"+(i+1)+"*"+coef[i] + "x^" + (i+1) + " : " + (1./(i+1) * coef[i])); + } + } + if(constI == null) + integ.constI = new double[]{1}; + else { +// System.out.println(degreeConstI()+2); + integ.constI = new double[degreeConstI()+1+1]; // degree +1 and +1 for beginning null in array +// for(int i = degreeConstI(); i >= 0; i--) { +// integ.constI[i+1] = 1./(i+1) * constI[i]; +// System.out.println(integ.constI[i]); +// System.out.println(integ.constI[i+1]); +// } + for(int i = 0; i <= degreeConstI(); i++) { + integ.constI[i+1] = 1./(i+1) * constI[i]; + } + integ.constI[0] = 1; +// System.out.println(integ.constI[0]); + } +// System.out.println(integ.constI[0]); + return integ; + } + + /** + * integrate this polynomial and return the value + * @param left + * @param right + * @return the value + */ + public double integrate(double left, double right){ + double result = 0; + for(int i=0; i<=degree(); i++){ + result += coef[i]/(i+1)*(Math.pow(right, (i+1)) - Math.pow(left, (i+1))); + } + return result; + } + + /** + * Include one condition for the integration constant of the constant term or + * for a polynomial with only one integration constant. + * @param position or value of the independent variable + * @param value of the condition + * @return the polynomial with the included condition + * @throws IllegalDimensionException + */ + public Polynomial condition(double position, double value) throws IllegalDimensionException{ + Polynomial result = new Polynomial(this); + // TODO: verify, seems good + if(degreeConstI() == 0 || (position == 0 && constI[0]>0) || countI()==1){ + int ic = 0; // the constant term case + if(countI()==1) ic = new Vector(constI).indexOfMax(); // the only one integration constant case + result.constI[ic] = 0; // delete the integration constant + if(result.getCoef(ic)>0) // if there is already a value + result.setCoef(ic,0); // integrate it to C (set the value to null / delete it) + result.setCoef((value-result.evaluate(position)/Math.pow(position,ic)), ic); + } else throw new IllegalDimensionException("Illegal integration constant dimensions."); + return result; + } + + /** + * Include exactly the same number of conditions as integrations constants. + * @param conditions triple of + * number of derivatives, + * the position or value of the independent variable and + * the value of the condition. + * @return the polynomial with the included conditions + * @throws IllegalDimensionException + * @throws SingularException + * @throws NoSquareException + */ + public Polynomial condition(double[][] conditions) throws IllegalDimensionException, NoSquareException, SingularException{ + int countI = countI(); + int countC = conditions.length; + if(countI != countC) + throw new IllegalDimensionException("The number of integration constants and the number of conditions must be the same."); + // TODO: verify, seems good + Matrix A = new Matrix(countI,countI); + Polynomial P = new Polynomial(this); + Vector b = new Vector(countI); + int i,j,k; + for(i=0;i0) P = this.derive(k); +// System.out.println(P); + if(P.constI!=null) + for(j=0;j0) + A.set(i, countI-1-k++, P.constI[j]*Math.pow(conditions[i][1],j)); +// System.out.println(conditions[i][2] + " - " + P.evaluate(conditions[i][1])); + b.set(i, conditions[i][2] - P.evaluate(conditions[i][1])); + } + Vector x = A.inverse().times(b); +// System.out.println(A); System.out.println(b); System.out.println(x); + + Polynomial result = new Polynomial(this); +// System.out.println(result.coef.length + " " + result.constI.length); + j = 0; + for(i=0;i2 + bx + c (where a != 0) + * will be replaced by x2 + px + q + * with p = b/a and q = c/a + * @return monic polynomial + */ + public Polynomial monicForm(){ + if(coef.length>1 && coef[coef.length-1]!=1){ + int i; + double[] monic = new double[coef.length]; + for(i=0; i + * Computation of cubic polynomials with Cardano's method if D>=0 for real roots, + * as a reminder:
+ * monic form x3+ax2+bx+c=0, substitute x = z -a/3
+ * z3+pz+q=0 with p=b-a2/3, q=2a3/27-a*b/3+c
+ * u=cbrt(-q/2+sqrt(D)), v=cbrt(-q/2-sqrt(D))
+ * e1=-1/2+1/2isqrt(3), e2=e12=-1/2-1/2isqrt(3)
+ * z1=u+v, z2=ue1+ve2, z3=ue2+ve1 + * @return the real roots/zeros of a function. + * @throws ComplexException + * @throws IllegalDimensionException + * @see Polynomial#roots() + */ + public Vector rootsRe() throws ComplexException, IllegalDimensionException{ + Vector roots = null; + VectorComplex cr = roots(); + switch(degree()){ + case 1: + roots = new Vector(cr.getRe(0)); + break; + case 2: + if(isC()) throw new ComplexException(); + roots = new Vector(cr.getRe(0),cr.getRe(1)); + break; + case 3: + if(!isC()) roots = new Vector(cr.getRe(0),cr.getRe(1),cr.getRe(2)); + else roots = new Vector(cr.getRe(0)); // x2, x3 complex + break; + default: + break; + } + if(degree() > 0) roots.sort(); + return roots; + } + + /** + * Calculate complex roots of a real function.
+ * Computation of cubic polynomials with Cardano's method if D>=0 for real roots, + * as a reminder:
+ * monic form x3+ax2+bx+c=0, substitute x = z -a/3
+ * z3+pz+q=0 with p=b-a2/3, q=2a3/27-a*b/3+c
+ * u=cbrt(-q/2+sqrt(D)), v=cbrt(-q/2-sqrt(D))
+ * e1=-1/2+1/2isqrt(3), e2=e12=-1/2-1/2isqrt(3)
+ * z1=u+v, z2=ue1+ve2, z3=ue2+ve1 + * @return the real roots/zeros of a function. + * @throws IllegalDimensionException + * @throws ComplexException + */ + public VectorComplex roots() throws IllegalDimensionException{ + if(degree()>3) throw new IllegalDimensionException(); + VectorComplex roots = null; + double a,b,c,D; + NumberComplex x1, x2 = null, x3; + switch(degree()){ + case 1: + roots = new VectorComplex(-coef[0]/coef[1]); + break; + case 2: + a = coef[2]; + b = coef[1]; + c = coef[0]; + D = b*b - 4*a*c; + if(D>=0){ + x1 = new NumberComplex((-b - Math.sqrt(D)) / (2*a)); + x2 = new NumberComplex((-b + Math.sqrt(D)) / (2*a)); + } else { + x1 = new NumberComplex(-b/(2*a),-Math.sqrt(-D)/(2*a)); + x2 = new NumberComplex(-b/(2*a), Math.sqrt(-D)/(2*a)); + } + roots = new VectorComplex(x1,x2); + break; + // TODO: validate !!! + case 3: + Polynomial monic = monicForm(); + a = monic.coef[2]; + b = monic.coef[1]; + c = monic.coef[0]; +// System.out.println(this); +// System.out.println(a+" "+b+" "+c); + double p = b - a*a/3; + double q = 2*a*a*a/27 - a*b/3 + c; +// System.out.println(p+" "+q); + D = q*q/4 + p*p*p/27; +// System.out.println("Polynomial:rootRe:D = " + D); + NumberComplex z1 = null, z2 = null, z3 = null; + if(p==0 && q==0){ + z1 = new NumberComplex(); + z2 = new NumberComplex(); + z3 = new NumberComplex(); + } + else if(D==0){ // verified! + double u = Math.cbrt(-q/2); // Math.cbrt(-q/2 +Math.sqrt(D)) + double v = Math.cbrt(-q/2); // Math.cbrt(-q/2 -Math.sqrt(D)) + z1 = new NumberComplex(u+v); + z2 = new NumberComplex(-(u+v)/2); // -(u+v)/2 -((u-v)/2)*Math.sqrt(3)*i + z3 = new NumberComplex(-(u+v)/2); // -(u+v)/2 +((u-v)/2)*Math.sqrt(3)*i + } + else if(D<0){ + double roh = Math.sqrt(-p*p*p/27); + double phi = Math.acos(-q/(2*roh)); + z1 = new NumberComplex(2*Math.cbrt(roh)*Math.cos(phi/3)); + z2 = new NumberComplex(2*Math.cbrt(roh)*Math.cos((phi+2*Math.PI)/3)); + z3 = new NumberComplex(2*Math.cbrt(roh)*Math.cos((phi+4*Math.PI)/3)); + } + else if(D>0){ // verified! + double u = Math.cbrt(-q/2 + Math.sqrt(D)); + double v = Math.cbrt(-q/2 - Math.sqrt(D)); + z1 = new NumberComplex(u + v); + z2 = new NumberComplex(-(u+v)/2, -((u-v)/2)*Math.sqrt(3)); + z3 = new NumberComplex(-(u+v)/2, +((u-v)/2)*Math.sqrt(3)); + } + x1 = z1.minus(a/3); + x2 = z2.minus(a/3); + x3 = z3.minus(a/3); + roots = new VectorComplex(x1,x2,x3); + break; + default: + break; + } + return roots; + } + + /** + * Run Newton's method to finds the roots of one variable. + * Newton's method to find x* such that f(x*) = 0, starting at x + * @param x + * @return roots of one variable + */ + public double rootNewton(double x) { + while(Math.abs(evaluate(x) / derive(x)) > Maths.ε) + x = x - evaluate(x) / derive(x); + return x; + } + + /** + * Run Newton's method to finds the local min/max of + * a twice-differentiable function of one variable. + * Newton's method to find x* such that f'(x*) = 0, starting at x + * @param x + * @return the local min/max + */ + public double optimumNewton(double x) { + while(Math.abs(derive(x) / derive(2,x)) > Maths.ε) + x = x - derive(x) / derive(2,x); + return x; + } + + /** + * Solve polynomial or computes intersections with polynomial + * @param b polynomial + * @return intersections + * @throws ComplexException + * @throws IllegalDimensionException + */ + public Vector solve(Polynomial b) throws ComplexException, IllegalDimensionException{ + return this.minus(b).rootsRe(); + } + + /** + * Check if this polynomial is complex. + * @return boolean true or false + */ + public boolean isC(){ + double D; + switch(degree()){ + case 1: + return true; + case 2: + D = Math.pow(coef[1], 2) - 4*coef[0]*coef[2]; + if(D>=0) return false; + else return true; + case 3: + Polynomial monic = monicForm(); + double a = monic.coef[2], b = monic.coef[1], c = monic.coef[0]; + double p = b - a*a/3, q = 2*a*a*a/27 - a*b/3 + c; + D = q*q/4 + p*p*p/27; + if(p==0 && q==0 || D<=0) return false; + else return true; + default: + return true; + } + } + + /** + * Is this Polynomial object equal to b? + * if(this == b) + * @param b g(x) + * @return boolean f(x) =? g(x) + */ + @Override + public boolean equals(Object b) { + if (b == null) return false; + if (b.getClass() != this.getClass()) return false; + Polynomial c = ((Polynomial) b).monicForm(); + Polynomial a = monicForm(); + if(a.degree() != c.degree()) return false; + int i; + for(i=0; i<=degree(); i++){ + if(a.coef[i] != c.coef[i]) return false; + } + return true; + } + + /** + * Convert to string representation + */ + @Override + public String toString() { + int i; + String s = ""; + double value; + String[] zahlaufteilung; + int pre = 1, post = 0, main = 0; + for(i = degree(); i >= 0; i--) { + value = (coef[i] <= Maths.ε && coef[i] >= -Maths.ε) ? 0.0 : coef[i]; + zahlaufteilung = String.valueOf(value).split("[.]"); // Am Komma trennen + if(zahlaufteilung.length>1 && !zahlaufteilung[1].equals("0")) + main = zahlaufteilung[1].length(); // max. Nachstellen + if(main != 0) pre++; // false = +1 für Komma + post = main>3 ? 3 : main; // 3 = max Nachstellen in der Ausgabe; +// System.out.println(Arrays.toString(zahlaufteilung) + " " + (pre+post)+"."+post); + if (i != degree() && value == 0) continue; + if (i == degree() && value < 0) s = s + "-"; + else if (i == degree()) s = s + ""; + else if (value > 0) s = s + " + "; + else if (value < 0) s = s + " - "; + if (i != 0 && (value == 1 || value == -1)) s = s + ""; +// if (i != 0 && value == 1) s = s + ""; +// else if (i != 0 && value ==-1) s = s + "-"; +// if (i == degree() || value < 0) s = s + "-"; + else if (i == degree() && value >= 0 || value > 0) s = s + String.format("%"+(pre+post)+"."+post+"f", value); + else if (value < 0) s = s + String.format("%"+(pre+post)+"."+post+"f",-value); + if (i == 1) s = s + symbolic; + else if (i > 1) s = s + symbolic+"^" + i; + } + + if(constI != null){ + String symbolicI = "c"; + int j = 1; + for(i = constI.length-1; i >= 0; i--) { + if(constI[i] > 0){ + s = s + " + "; + if (constI[i] != 1) s = s + ( constI[i] ); + s = s + symbolicI + j++; + if (i == 1) s = s + symbolic; + else if (i > 1) s = s + symbolic+"^" + i; + } + } + } + return s; + } + + public String toStr(){ + StringBuffer output = new StringBuffer(); + for(int i = degree(); i >= 0; i--) output.append(coef[i]+" "); + return output.toString(); + } + + public void println(){ + System.out.println(toString()); + } + + public void plot(){ + Vector X = Vector.fill(0, 0.001, 1d); + Vector Y = new Vector(evaluate(X.getArray())).times(1); + Draw draw = new Draw(); + draw.animate(X, Y); + } + + /** + * test client + * @param args + */ + public static void main(String[] args) { + Polynomial zero = new Polynomial(0, 0); + Polynomial con0 = new Polynomial(5, 0); + + Polynomial p1 = new Polynomial(4, 3); + Polynomial p2 = new Polynomial(3, 2); + Polynomial p3 = new Polynomial(2, 0); + Polynomial p4 = new Polynomial(1, 1); + + // all 4 lines do the same thing + Polynomial p = p1.plus(p2).plus(p3).plus(p4); // 4x^3 + 3x^2 + x + 2 + p = new Polynomial(new double[]{4,3,1,2}); + p = new Polynomial(" 4 3 1 2"); + p = new Polynomial(4,3,1,2); + + Polynomial q1 = new Polynomial(3, 2); + Polynomial q2 = new Polynomial(5, 0); + Polynomial q = q1.plus(q2); // 3x^2 + 5 + + Polynomial r = new Polynomial("4 +8 1 -1"); + + System.out.println("zero(x) = " + zero); + System.out.println("p(x) = " + p); + System.out.println("p(0) = " + p.evaluate(0)); + System.out.println("q(x) = " + q); + System.out.println("p(x) + q(x) = " + p.plus(q)); + System.out.println("p(x) * q(x) = " + p.times(q)); + System.out.println("p(q(x)) = " + p.compose(q)); + System.out.println("0 - p(x) = " + zero.minus(p)); + System.out.println("p(3) = " + p.evaluate(3)); + System.out.println("p(x) = " + p); + System.out.println("p'(x) = " + p.derive()); + System.out.println("p''(x) = " + p.derive().derive()); + System.out.println("monic p''(x)= " + p.derive().derive().monicForm()); + System.out.println("-1 * p''(x) = " + p.derive().derive().monicForm().times(-1)); + System.out.println("q(x) = " + q); + System.out.println("Q(1,5) = " + q.integrate(1, 5)); + System.out.println("Q = " + q.integrate()); + System.out.println("int(Q) = " + q.integrate().integrate()); + try { + System.out.println("int(Q(0))=5 = " + q.integrate().integrate().condition(0, 5)); + System.out.println("int(int(Q)) = " + q.integrate().integrate().condition(0, 5).integrate()); + } catch (IllegalDimensionException e) { + e.printStackTrace(); + } + System.out.println("deg(q) = " + q.degree()); + System.out.println("a = " + con0); + System.out.println("a' = " + con0.derive()); + System.out.println("p(x) = " + p); + try { + System.out.println("roots(p) = \n" + p.rootsRe()); + System.out.println("r(x) = " + r); + System.out.println("roots(r) = \n" + r.rootsRe()); + } catch (ComplexException | IllegalDimensionException e1) { + // TODO Auto-generated catch block + e1.printStackTrace(); + } + System.out.println("p =? q = " + p.equals(q)); + Polynomial w = new Polynomial(-1, 0).integrate().integrate().integrate().integrate(); + System.out.println("w = " + w); + System.out.println("w'''' = " + w.derive().derive().derive().derive()); + System.out.println("w'''' = " + w.derive(4)); + try { + System.out.println("w = " + w.condition(new double[][]{ + {0,0,0},{1,0,0},{2,1,0},{3,1,0} + })); + } catch (IllegalDimensionException | NoSquareException | SingularException e) { + e.printStackTrace(); + } +// Vector X = Vector.fill(0, 0.001, 1d); +// Vector Y = new Vector(w.evaluate(X.get())).times(1); +// Draw draw = new Draw(); +// draw.animate(X, Y); + + System.out.println(); + Polynomial T9 = Polynomial.chebyshevT(9); + System.out.println(T9); + Polynomial T9T = new Polynomial(256,0,-576,0,432,0,-120,0,9,0); + System.out.println("Chebyshev T(9) test: "+T9.equals(T9T)); + Polynomial U9 = Polynomial.chebyshevU(9); + System.out.println(U9); + Polynomial U9T = new Polynomial(512,0,-1024,0,672,0,-160,0,10,0); + System.out.println("Chebyshev U(9) test: "+U9.equals(U9T)); + + Polynomial H9 = Polynomial.hermite(9); + System.out.println(H9); + Polynomial H9T = new Polynomial(512,0,-9216,0,48384,0,-80640,0,30240,0); + System.out.println("Hermite H(9) test: "+H9.equals(H9T)); + Polynomial H9d = Polynomial.hermiteD(9); + System.out.println(H9d); + Polynomial H9dT = new Polynomial(4608,0,-64512,0,241920,0,-241920,0,30240); + System.out.println("Hermite H'(9) test: "+H9d.equals(H9dT)); + Polynomial H9d2 = Polynomial.hermiteD(9,2); + System.out.println(H9d2); + Polynomial H9d2T = new Polynomial(36864,0,-387072,0,967680,0,-483840,0); + System.out.println("Hermite H''(9) test: "+H9d2.equals(H9d2T)); + + Polynomial He9 = Polynomial.hermiteE(9); + System.out.println(He9); + Polynomial He9T = new Polynomial(1,0,-36,0,378,0,-1260,0,945,0); + System.out.println("Hermite He(9) test: "+He9.equals(He9T)); + Polynomial He9d = Polynomial.hermiteED(9); + System.out.println(He9d); + Polynomial He9dT = new Polynomial(9,0,-252,0,1890,0,-3780,0,945); + System.out.println("Hermite He'(9) test: "+He9d.equals(He9dT)); + Polynomial He9d2 = Polynomial.hermiteED(9,2); + System.out.println(He9d2); + Polynomial He9d2T = new Polynomial(72,0,-1512,0,7560,0,-7560,0); + System.out.println("Hermite He''(9) test: "+He9d2.equals(He9d2T)); + } + +} \ No newline at end of file diff --git a/src/math/equation/PolynomialComplex.java b/src/math/equation/PolynomialComplex.java new file mode 100644 index 0000000..efbb3e0 --- /dev/null +++ b/src/math/equation/PolynomialComplex.java @@ -0,0 +1,114 @@ +package math.equation; + +import math.number.NumberComplex; + + +public class PolynomialComplex { + private Polynomial re; + private Polynomial im; + + public PolynomialComplex(){ + setRe(new Polynomial()); + im = new Polynomial(); + } + + public PolynomialComplex(String cpol){ + this(new Polynomial(cpol)); + } + + public PolynomialComplex(Polynomial c){ + int i; + for(i=0; i<=c.degree(); i++){ + // TODO: + } + + int size = c.degree()+1; + double[] coefRe = new double[size]; + double[] coefIm = new double[size]; + for(i=0; i q.length() ? p.length() : q.length(); + String sep = ""; + for(int i=0; i q.length() ? p.length() : q.length(); + String sep = ""; + for(int i=0; ii∈ℂ + * @author Daniel + * + */ +public class VectorComplex extends WObject{ + /** + * ai0 = real ℜ, + * ai1 = imaginary ℑ + */ + private double[][] vector; + + public VectorComplex(){ + } + + /** + * Create zero vector + * @param n rows + */ + public VectorComplex(int n){ + vector = new double[n][2]; + int i; + for(i=0; i0 ? new NumberComplex(vector[0][0],vector[0][1]) : null; + } + + /** + * get y (second element) + * @return y value + */ + public NumberComplex y(){ + return n()>1 ? new NumberComplex(vector[1][0],vector[1][1]) : null; + } + + /** + * get z (third element) + * @return z value + */ + public NumberComplex z(){ + return n()>2 ? new NumberComplex(vector[2][0],vector[2][1]) : null; + } + + /** + * Create vector with zeros + * @param n size + * @return zero vector + */ + public static VectorComplex zeros(int n){ + return new VectorComplex(n); + } + + /** + * Create vector with ones + * @param n size + * @return one vector + */ + public static VectorComplex ones(int n){ + return fill(n, 1); + } + + /** + * Create vector with given complex number + * @param n size + * @param s real number + * @return complex vector filled with real numbers + */ + public static VectorComplex fill(int n, double s){ + return fill(n, s, 0); + } + + /** + * Create vector with given real and imaginary part + * @param n size + * @param r real + * @param i imaginary + * @return vector filled with complex numbers + */ + public static VectorComplex fill(int n, double r, double i){ + VectorComplex a = new VectorComplex(n); + int j; + for(j=0; j1 ? 1 : deltaRe<-1 ? -1 : 0, + deltaIm>1 ? 1 : deltaIm<-1 ? -1 : 0); + return fill(start, increment, end); + } + + /** + * Create vector and fill it from a start value and increments it to an end value + * @param start + * @param increment + * @param end + * @return filled vector [start,start+increment,start+2*increment,...,end] + */ + public static VectorComplex fill(NumberComplex start, NumberComplex increment, NumberComplex end){ + double deltaRe = end.ℜ()-start.ℜ(); + double deltaIm = end.ℑ()-start.ℑ(); + int nRe = increment.ℜ()!=0?(int)(deltaRe/increment.ℜ())+1:0; + int nIm = increment.ℑ()!=0?(int)(deltaIm/increment.ℑ())+1:0; + int n = nRe > nIm ? nRe : nIm; + if(n==0||(deltaRe>=0&&increment.ℜ()<0&&deltaIm>=0&&increment.ℑ()<0)|| + (deltaRe<0&&increment.ℜ()>=0&&deltaIm<0&&increment.ℑ()>=0)|| + (increment.ℜ()==0&&increment.ℑ()==0)) return new VectorComplex(start); + VectorComplex a = new VectorComplex(n); + int i; + for(i=0; i 0){ + result = false; + break; + } + if(vector[i][1] > 0){ + result = false; + break; + } + } + return result; + } + + @Override + public String toString(){ + int i; + String output = ""; + if(vector!=null){ + for(i=0; i-1E-10) ? null : vector[i]; + output += String.format("% 10.4g % 10.4gi \n", vector[i][0], vector[i][1]); + } + } + return output; + } + + /** + * Transpose (horizontal) + * @return string of vector in transposed form + */ + public String toStringT(){ + int i; + String output = ""; + if(vector!=null){ + for(i=0; i-1E-10) ? null : vector[i]; + output += String.format("% 10.4g % 10.4gi ", vector[i][0], vector[i][1]); + } + output += "\n"; + } + return output; + } + + /** + * Transpose (horizontal) + */ + public void toConsoleT(){ + System.out.println(toStringT()); + } + + /** + * test client + * @param args + */ + public static void main(String[] args) { + VectorComplex a = new VectorComplex(new NumberComplex(5), new NumberComplex(1), new NumberComplex(1)); + System.out.println(a.toString()); + VectorComplex.fill(new NumberComplex(2.), new NumberComplex(2.5)).println(); + VectorComplex.fill(new NumberComplex(2.), new NumberComplex(8)).println(); + VectorComplex.fill(new NumberComplex(8.), new NumberComplex(2)).println(); + VectorComplex.fill(new NumberComplex(2.), new NumberComplex(2), new NumberComplex(8)).println(); + VectorComplex.fill(new NumberComplex(8.), new NumberComplex(-2), new NumberComplex(2)).toConsoleT(); + System.out.println("a isNull? : " + a.isZero()); + System.out.println("a-a isNull? : " + a.minus(a).isZero()); + } +} diff --git a/src/math/number/Number.java b/src/math/number/Number.java new file mode 100644 index 0000000..cbd31af --- /dev/null +++ b/src/math/number/Number.java @@ -0,0 +1,12 @@ +package math.number; + +import thisandthat.WObject; + +public class Number extends WObject { + + public static void main(String[] args) { + // TODO Auto-generated method stub + + } + +} diff --git a/src/math/number/NumberComplex.java b/src/math/number/NumberComplex.java new file mode 100644 index 0000000..c757b5a --- /dev/null +++ b/src/math/number/NumberComplex.java @@ -0,0 +1,377 @@ +package math.number; + +import math.Maths; + + +/** + * Data type for complex number ∈ℂ which defines complex arithmetic and + * mathematical functions. + * @author Daniel Weschke + */ +public final class NumberComplex extends Number{ + + /** + * real part ℜ + */ + private double ℜ; + + /** + * imaginary part ℑ + */ + private double ℑ; + + /** + * Constructs the complex number z = a + bi, a,b=0 + */ + public NumberComplex(){ + } + + /** + * Constructs the complex number z = a + bi, b=0 + * @param a real part + */ + public NumberComplex(double a){ + ℜ = a; + } + + /** + * Constructs the complex number z = ℜ(z) + ℑ(z)i + * @param ℜ real part + * @param ℑ imaginary part + */ + public NumberComplex(double ℜ, double ℑ){ + this.ℜ = ℜ; + this.ℑ = ℑ; + } + + /** + * Constructs the complex number z = a + bi + * @param z complex number + */ + public NumberComplex(NumberComplex z){ + this(z.ℜ(), z.ℑ()); + } + + /** + * Constructs the complex number z = a + bi + * @param z complex number + * @return complex number + */ + public NumberComplex assign(NumberComplex z){ + ℜ = z.ℜ(); + ℑ = z.ℑ(); + return this; + } + + /** + * Real part of this complex number + * (the x-coordinate in rectangular coordinates). + * @return Re(z) where z is this complex number + */ + public double ℜ(){ + return ℜ; + } + + /** + * Imaginary part of this complex number + * (the y-coordinate in rectangular coordinates). + * @return ℑ(z) where z is this complex number + */ + public double ℑ(){ + return ℑ; + } + + /** + * Addition with real number (doesn't change this complex number).
+ * (a+bi) + c = (a+c)+bi. + * @param c real number + * @return (a+bi) + c + * @see #plus(NumberComplex) + * @see #plusAssign(NumberComplex) + */ + public NumberComplex plus(double c){ + return new NumberComplex(ℜ()+c,ℑ()); + } + + /** + * Addition with complex number (doesn't change this complex number).
+ * (a+bi) + (c+di) = (a+c)+(b+d)i. + * @param z complex number + * @return (a+bi) + (c+di) + * @see #plus(double) + * @see #plusAssign(NumberComplex) + */ + public NumberComplex plus(NumberComplex z){ + return new NumberComplex(ℜ()+z.ℜ(),ℑ()+z.ℑ()); + } + + /** + * Addition with complex number (does change this complex number).
+ * (a+bi) + (c+di) = (a+c)+(b+d)i. + * @param z complex number + * @return y = y + z + * @see #plus(double) + * @see #plus(NumberComplex) + */ + public NumberComplex plusAssign(NumberComplex z){ + ℜ += z.ℜ(); + ℑ += z.ℑ(); + return this; + } + + /** + * Subtraction of Complex numbers (doesn't change this Complex number).
+ * (a+bi) - c = (a-c)+bi. + * @param c real number + * @return (a+bi) - c + * @see #minus(NumberComplex) + */ + public NumberComplex minus(double c){ + return plus(-c); + } + + /** + * Subtraction of Complex numbers (doesn't change this Complex number).
+ * (a+bi) - (c+di) = (a-c)+(b-d)i. + * @param z complex number + * @return (a+bi) - (c+di) + * @see NumberComplex#minus(double) + */ + public NumberComplex minus(NumberComplex z){ + return new NumberComplex(ℜ()-z.ℜ(),ℑ()-z.ℑ()); + } + + /** + * Scalar multiplication (doesn't change this Complex number). + * @param s scalar + * @return s(a+bi) + * @see #times(NumberComplex) + */ + public NumberComplex times(double s){ + return new NumberComplex(ℜ()*s,ℑ()*s); + } + + /** + * Complex multiplication (doesn't change this Complex number). + * @param z complex number + * @return (a+bi)(c+di) = ac-bd + (ad+bc)i + * @see #times(double) + */ + public NumberComplex times(NumberComplex z){ + // local r and i is needed because re and im depends on re and im + double r = ℜ()*z.ℜ() - ℑ()*z.ℑ(); + double i = ℜ()*z.ℑ() + ℑ()*z.ℜ(); + return new NumberComplex(r, i); + } + + /** + * Division of Complex numbers (doesn't change this Complex number).
+ * (a+bi)/(c+di) = (ac+bd)/(cc+dd) + (bc-ad)/(cc+dd) i + * @param z complex number + * @return (a+bi)/(c+di) = (ac+bd)/(cc+dd) + (bc-ad)/(cc+dd) i + */ + public NumberComplex over(NumberComplex z){ + return times(z.reciprocal()); + } + + /** + * @return a new complex object whose value is the reciprocal of this + */ + public NumberComplex reciprocal() { + double scale = ℜ()*ℜ() + ℑ()*ℑ(); + return new NumberComplex(ℜ() / scale, -ℑ() / scale); + } + + /** + * Modulus or absolute value or length of this complex number + * (the distance from the origin in polar coordinates). + * @return |z| where z is this Complex number. + */ + public double mod(){ + return Maths.hypot(ℜ(), ℑ()); + } + + /** + * Argument or phase of this Complex number + * (the angle in radians with the real/x-axis in polar coordinates). + * @return arg(z) in radians where z is this Complex number + */ + public double arg(){ + return Maths.atan(ℑ(), ℜ()); + } + + /** + * Argument of this complex number + * (the angle ∠ in radians with the real/x-axis in polar coordinates). + * @return arg(z) in deg where z is this complex number + */ + public double argd(){ + return Maths.atand(ℑ(), ℜ()); + } + + /** + * Complex conjugate of this complex number + * (the conjugate of x+i*y is x-i*y). + * @return z-bar where z is this complex number. + */ + public NumberComplex conj() { + return new NumberComplex(ℜ(),-ℑ()); + } + + /** + * Complex square root (doesn't change this complex number). + * Computes the principal branch of the square root, which + * is the value with 0 <= arg < pi. + * @return sqrt(z) where z is this Complex number. + */ + public NumberComplex sqrt(){ + double r=Math.sqrt(this.mod()); + double ϑ=this.arg()/2; + return new NumberComplex(r*Math.cos(ϑ),r*Math.sin(ϑ)); + } + + /** + * Complex exponential (doesn't change this Complex number). + * @return exp(z) where z is this Complex number. + */ + public NumberComplex exp() { + return new NumberComplex(Math.exp(ℜ())*Math.cos(ℑ()), + Math.exp(ℜ())*Math.sin(ℑ())); + } + + /** + * Principal branch of the complex logarithm of this complex number + * (doesn't change this complex number). + * The principal branch is the branch with -π < arg <= π. + * @return log(z) where z is this Complex number. + */ + public NumberComplex log() { + return new NumberComplex(Math.log(this.mod()),this.arg()); + } + + /** + * Cosine of this complex number (doesn't change this complex number).
+ * cos(z) = (eiz+e-iz)/2. + * @return cos(z) where z is this Complex number. + */ + public NumberComplex cos(){ + return new NumberComplex(cosh(ℑ())*Math.cos(ℜ()),-sinh(ℑ())*Math.sin(ℜ())); + } + + /** + * Sine of this complex number (doesn't change this complex number).
+ * sin(z) = (eiz-e-iz)/(2i). + * @return sin(z) where z is this Complex number. + */ + public NumberComplex sin() { + return new NumberComplex(cosh(ℑ())*Math.sin(ℜ()),sinh(ℑ())*Math.cos(ℜ())); + } + + /** + * Tangent of this complex number (doesn't change this complex number).
+ * tan(z) = sin(z)/cos(z). + @return tan(z) where z is this Complex number. + */ + public NumberComplex tan() { + return sin().over(cos()); + } + + /** + * Real cosh function (used to compute complex trig functions) + * @param ϑ argument + * @return (eϑ+e-ϑ)/2 + */ + private double cosh(double ϑ){ + return (Math.exp(ϑ)+Math.exp(-ϑ))/2; + } + + /** + * Hyperbolic cosine of this complex number (doesn't change this complex number).
+ * cosh(z) = (ez + e-z)/2. + * @return cosh(z) where z is this Complex number. + */ + public NumberComplex cosh() { + return new NumberComplex(cosh(ℜ())*Math.cos(ℑ()),sinh(ℜ())*Math.sin(ℑ())); + } + + /** + * Real sinh function (used to compute complex trig functions) + * @param ϑ argumant + * @return (eϑ-e-ϑ)/2 + */ + private double sinh(double ϑ){ + return (Math.exp(ϑ)-Math.exp(-ϑ))/2; + } + + /** + * Hyperbolic sine of this complex number (doesn't change this complex number).
+ * sinh(z) = (ez-e-z)/2. + * @return sinh(z) where z is this Complex number. + */ + public NumberComplex sinh(){ + return new NumberComplex(sinh(ℜ())*Math.cos(ℑ()),cosh(ℜ())*Math.sin(ℑ())); + } + + /** + * Negative of this complex number (chs stands for change sign). + * This produces a new Complex number and doesn't change this Complex number.
+ * -(x+i*y) = -x-i*y. + * @return -z where z is this Complex number. + */ + public NumberComplex chs() { + return new NumberComplex(-ℜ(),-ℑ()); + } + + /** + * A string representation of the complex object. + */ + @Override + public String toString(){ + if(ℑ() == 0) return ""+ℜ(); + if(ℜ() == 0) return ℑ()+"i"; + if(ℑ() > 0) return ""+ℜ()+" + "+ℑ()+"i"; + if(ℑ() < 0) return ""+ℜ()+" - "+(-ℑ())+"i"; + else return ℜ()+" + i*"+ℑ(); // shouldn't get here (unless Inf or NaN) + } + + public void printRe(){ + System.out.println("ℜ(z) = "+ℜ()); + } + + public void printIm(){ + System.out.println("ℑ(z) = "+ℑ()); + } + + /** + * @param args + */ + public static void main(String[] args) { + NumberComplex cn1 = new NumberComplex(4,1); + System.out.println("a = "+cn1); + NumberComplex cn2 = new NumberComplex(4,-1); + System.out.println("b = "+cn2); + NumberComplex cn = cn1.times(cn2); + System.out.println("a * b = "+cn); + System.out.println("mod(): "+(new NumberComplex(4,3).mod()==5)); + System.out.println(cn.argd()); + + + NumberComplex a = new NumberComplex(5.0, 6.0); + NumberComplex b = new NumberComplex(-3.0, 4.0); + + System.out.println("a = " + a); + System.out.println("b = " + b); + System.out.println("ℜ(a) = " + a.ℜ()); + System.out.println("ℑ(a) = " + a.ℑ()); + System.out.println("b + a = " + b.plus(a)); + System.out.println("a - b = " + a.minus(b)); + System.out.println("a * b = " + a.times(b)); + System.out.println("b * a = " + b.times(a)); + System.out.println("a / b = " + a.over(b)); + System.out.println("(a / b) * b = " + a.over(b).times(b)); + System.out.println("conj(a) = " + a.conj()); + System.out.println("|a| = " + a.mod()); + System.out.println("tan(a) = " + a.tan()); + } + +} diff --git a/src/physics/Nyquist.java b/src/physics/Nyquist.java new file mode 100644 index 0000000..65b42c8 --- /dev/null +++ b/src/physics/Nyquist.java @@ -0,0 +1,154 @@ +package physics; + +import java.awt.event.ActionEvent; +import java.awt.event.ActionListener; +import java.util.ArrayList; + +import javax.swing.JTextField; + +import math.Maths; +import math.equation.RationalPolynomial; +import math.equation.RationalPolynomialComplex; +import math.matrix.Matrix; +import math.matrix.Vector; +import stdlib.StdDraw; +import awt.Draw; + +public class Nyquist{ + Draw nyquist = new Draw(); + static JTextField tf = new JTextField(); + static Thread plot; + static Vector X; + static Vector Y; + + public Matrix getNyquist() { + RationalPolynomial PT11 = new RationalPolynomial("2; 4 1"); + RationalPolynomial PT12 = new RationalPolynomial("1; 6 1"); +// RationalFunction I = new RationalFunction("5; 24 0"); + RationalPolynomial sys = PT11.times(PT12).times(PT12).times(PT12); + tf.setText(sys.a.toStr()+"; "+sys.b.toStr()); + RationalPolynomialComplex crf = new RationalPolynomialComplex(sys); +// doule v = crf.getDivident().getRe().getCoef(0); + PT11.setSymbolic("s"); + PT12.setSymbolic("s"); + sys.setSymbolic("s"); + crf.setSymbolic("s"); + return getNyquist(crf); + } + + public Matrix getNyquist(RationalPolynomial sys) { + tf.setText(sys.a.toStr()+"; "+sys.b.toStr()); + RationalPolynomialComplex crf = new RationalPolynomialComplex(sys); + return getNyquist(crf); + } + + public Matrix getNyquist(RationalPolynomialComplex crf){ + ArrayList X = new ArrayList(); + ArrayList Y = new ArrayList(); + double omega, abs, angle; + double x=0, y=0, t=0; + double epsilon = 1; + double directionX, directionY; + while(epsilon > 0.01){ + omega = t; + abs = crf.mod(omega); + angle = crf.argd(omega); + directionX = + abs * Maths.cosd(angle) - x; + directionY = - abs * Maths.sind(angle) - y; + + // Determine the pointers location + x += directionX; + y += directionY; + t += 0.001; + + X.add(x); + Y.add(y); + epsilon = Maths.hypot(x, y); + } + + // get array + Object xa[] = X.toArray(); + Object ya[] = Y.toArray(); + int imax = xa.length; + double[] xd = new double[imax]; + double[] yd = new double[imax]; + // the array + for(int i=0; i doppel thead? + Draw.clear(); + nyquist.show(1); + Draw.setProcess(0); + RationalPolynomial sys = new RationalPolynomial(tf.getText().toString()); + setNyquist(sys); + animate(); + } + }); + } + + public void animate(){ + if(X==null) demo(); + nyquist.markX(-1); + nyquist.animate(); + } + + /** + * 2.0 ; 264.0 348.0 80.0 22.0 1.0 + * 10.0 2.0 1.0 ; 864.0 648.0 180.0 22.0 + * 2.0 ; 64.0 1.0 1.0 + * @param args + */ + public static void main(String[] args) { + Nyquist ny = new Nyquist(); +// ny.setNyquist(new RationalPolynomial("10.0 2.0 1.0 ; 864.0 648.0 180.0 22.0")); +// ny.setNyquist(new RationalPolynomial("10.0 2.0 1.0 ; 864.0 648.0 180.0 22.0")); + ny.animate(); + ny.nyquist.cross(); + System.out.println(); + } +}