Petite librairie d'analyse numérique

Version imprimable

Je vous présente ci-dessous une petite librairie d'analyse numérique, comprenant les domaines suivants :

- Recherche de racines
- Résolution de sysètmes linéaires
- Interpolation
- Intégration

Le code se veut souple et adaptable, je poste d'abord les sources, les examples d'utilisations suivront, au cas où la javadoc ne suffirait pas.

Voici d'abord quelques classes d'ordre général utilisées en analyse et dans les domaines présentés ci-dessus :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 
/*
 * Function.java
 *
 * Created on 28 novembre 2006, 10:25
 *
 */
 
package org.umh.math.analysis;
 
/**
 * A function is represented by its evaluation method.
 * @author Absil Romain
 */
public interface Function
{
    /**
     * Evaluate the function on the given argument.
     * @return the evaluation of the function on the given argument.
     * @param x the number you want the function to be evaluated.
     */
    public double eval(double x);
}

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
 
/*
 * Function2D.java
 *
 * Created on 16 janvier 2007, 13:42
 *
 */
 
package org.umh.math.analysis;
 
import java.awt.geom.Point2D;
 
/**
 * This class models functions from R to R² wich are represented by 
 * their evaluation.
 * @author Absil Romain
 */
public class Function2D
{
    private Function fy;
    private Function fx;
 
    /**
     * Construcs a new Function2D with the functions defining its evaluation.
     * @param fx the function giving the x coordinate of the evaluations.
     * @param fy the function giving the y coordinate of the evaluations.
     **/
    public Function2D(Function fx, Function fy)
    {
        this.fx = fx;
        this.fy = fy;
    }
 
    /**
     * Returns the function giving the y coordinate of the evaluations.
     * @return the function the y coordinate of the evaluations.
     **/
    public Function getFy()
    {
        return fy;
    }   
 
    /**
     * Sets the function giving the y coordinate of the evaluations.
     * @param fy the function to set.
     **/
    public void setFy(Function fy)
    {
        this.fy = fy;
    }
 
    /**
     * Returns the function giving the x coordinate of the evaluations.
     * @return the function giving the x coordinate of the evaluations.
     **/
    public Function getFx()
    {
        return fx;
    }
 
    /**
     * Sets the function giving the x coordinate of the evaluations.
     * @param fx the function to set.
     **/
    public void setFx(Function fx)
    {
        this.fx = fx;
    }
 
    /**
     * Evaluates the function to the given parameter.
     * @param t the parameter you want the function to be evaluated on.
     * @return the evaluation of the function to the given parameter.
     **/
    public Point2D.Double eval(double t)
    {
        return new Point2D.Double(fx.eval(t),fy.eval(t));
    }
}

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 
/*
 * Function_1_2.java
 *
 * Created on 22 mars 2007, 15:55
 *
 */
 
package org.umh.math.analysis;
 
import java.awt.geom.Point2D;
 
/**
 * This class models functions from R² to R.
 * @author Absil Romain
 */
public interface Function_2_1
{
    /**
     * Evaluates the function to the given parameters x and y.
     * @param x the first number you want the function to be evaluated on.
     * @param y the second number you want the function to be evaluated on.
     * @return the function evaluated to the point (x,y).
     */
    public double eval(double x, double y);
}

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
 
/*
 * Polynome.java
 *
 * Created on 22 mars 2007, 16:26
 *
 */
 
package org.umh.math.analysis;
 
/**
 * A polynome is a function on wich some other operations are defined.
 * @author Absil Romain
 */
public class Polynome implements Function
{
    private int deg;
    private double[] coffs;
 
    /**
     * Constructs a null polynome of degree n.
     * @param deg the degree of the polynome.
     **/
    public Polynome(int deg)
    {
        if(deg < 0)
            throw new IllegalArgumentException();
        this.deg = deg;
        this.coffs = new double[deg+1];
    }
 
    /**
     * Construcs a new polynome wich the coefficients are given in their
     * natural order.
     * <br>e.g. the call <code>new Polynome(1,2,3);</code> will creates the
     * polynome x² + 2x + 3.
     * @param coffs the coefficients of the polynome.
     **/
    public Polynome(double ... coffs)
    {
        this.deg = coffs.length - 1;
        this.coffs = new double[coffs.length];
        for (int i = 0; i < coffs.length; i++)
            this.coffs[i] = coffs[coffs.length - 1 - i];
        trim();
    }
 
    /**
     * Returns the degree of the polynome.
     * @return the degree of the polynome.
     **/
    public int getDeg()
    {
        return deg;
    }
 
    /**
     * Returns all the coefficients of the polynome. The coefficient indexed 
     * by the integer i is the coefficient of degree i.
     * @return all the coefficients of the polynome. The coefficient indexed 
     * by the integer i is the coefficient of degree i.
     **/
    public double[] getCoffs()
    {
        return coffs;
    }
 
    /**
     * Returns the coefficient of the given degree.
     * @param deg the degree of the coefficient you want to return.
     * @return the coefficient of the given degree.
     **/
    public double getCof(int deg)
    {
        return coffs[deg];
    }
 
    /**
     * Sets the coefficient of the given degree to the given value.
     * @param deg the degree of the coefficient yoy want to set.
     * @param value the coefficient to set.
     **/
    public void setCof(int deg, double value)
    {
        coffs[deg] = value;
    }
 
    /**
     * Sets all the coefficients of the current polynome by other coefficients.
     * @param coffs the coefficients to set.
     **/
    public void setCoffs(double[] coffs)
    {
        this.coffs = coffs;
    }
 
    /**
     * Returns the String representation of this polynome.
     * @return the String representation of this polynome.
     **/
    public String toString()
    {
        if(deg == 0 )
            return "" + coffs[0];
 
	String s = "";
	for(int i = deg ; i >= 0 ; i--)
	    {
		if(i != deg)
		    {
			if(getCof(i) > 0 && i != 1 && i != 0)
			    s += " + " + getCof(i) + " x^" + i;
			else if(getCof(i) < 0 && i != 3 && i !=  2 && i != 1 && i != 0)
			    s += " - " +  Math.abs(getCof(i)) + " x^" + i;
			else if(getCof(i) > 0 && i == 1)
			    s += " + " + getCof(i) + " x";
			else if(getCof(i) > 0 && i == 0)
			    s += " + " + getCof(i);
			else if(getCof(i) < 0 && i == 1)
			    s += " - " + Math.abs(getCof(i)) + " x";
			else if(getCof(i) < 0 && i == 0)
			    s += " - " + Math.abs(getCof(i));
		    }
		else if(i == deg)
		    {
			if(getCof(i) > 0 && i > 1)
			    s += getCof(i) + " x^" + i;
			else if(getCof(i) < 0 && i > 3)
			    s += Math.abs(getCof(i)) + " x^" + i;
			else if(getCof(i) > 0 && i == 1)
			    s += getCof(i) + " x";
			else if(getCof(i) > 0 && i == 0)
			    s += getCof(i);
			else if(getCof(i) < 0 && i == 1)
			    s += Math.abs(getCof(i)) + " x";
			else if(getCof(i) < 0 && i == 0)
			    s += Math.abs(getCof(i));
                        else if(getCof(i) == 0)
                            s += "";
		    }
	    }
	return s;
    }
 
    /**
     * Returns true if the given Object is equal to the current polynome,
     * returns false otherwise.
     * @return true if the given Object is equal to the current polynome,
     * returns false otherwise.
     * @param other the bject you want to know if it is equal to the current Polynome
     */
    public boolean equals(Object other)
    {
        if(!(other instanceof Polynome))
            return false;
        Polynome pol = (Polynome)other;
        pol.trim();
        if(pol.getDeg() != deg)
            return false;
        for (int i = 0; i < coffs.length; i++)
            if(coffs[i] != pol.getCof(i))
                return false;
        return true;
    }
 
    /**
     * Returns a copy of the current polynome.
     * @return a copy of the current polynome.
     **/
    public Polynome clone()
    {
        return new Polynome(coffs);
    }
 
    /**
     * Deletes all the first null coefficiets and update the degree of the 
     * polynome.
     **/        
    public void trim()
    {
        int i = deg;
        for(; i >= 0 ; i--)
        {
            if(coffs[i] != 0)
                break;            
        }
 
        if(i != deg)
        {
            double[] coffTmp = new double[i+1];
            for(int j = 0 ; j <= i ; j++)        
                coffTmp[j] = coffs[j];    
            deg = i;
            coffs = coffTmp;
        }        
    }
 
    /**
     * Returns the image of the polynome evaluated in x.
     * @param x the value you want to know the image.
     * @return the image of the polynome evaluated in x.
     **/
    public double eval(double x)
    {
	double img = 0;
        for (int i = deg; i >= 0; i--)
            img = coffs[i] + x * img;
        return img;
    }
 
    /**
     * Adds a polynome to the current polynome.
     * @param polynome the polynome to add.
     * @return the result of the addition.
     **/
    public Polynome add(Polynome polynome)
    {
	int newDeg = Math.max(deg,polynome.getDeg());
 
        //sets the degree of the 2 polynomes to add
	Polynome thisTmp = new Polynome(newDeg);        
	for(int i = 0 ; i <= deg ; i++)
	    thisTmp.setCof(i,getCof(i));
 
	Polynome polTmp = new Polynome(newDeg);
	for(int i = 0 ; i <= polynome.getDeg() ; i++)
	    polTmp.setCof(i,polynome.getCof(i));
 
        //adds the polynomes
	for(int i = 0 ; i <= newDeg ; i++)
	    thisTmp.setCof(i,thisTmp.getCof(i) + polTmp.getCof(i));
        thisTmp.trim();
        return thisTmp;
    }        
 
    /**
     * Returns the opposite of the current polynome for the additive law.
     * @return the opposite of the current polynome for additive law.
     */
    public Polynome getOpposite()
    {
        double[] newcofs = new double[coffs.length];
        for (int i = 0; i < newcofs.length; i++)
            newcofs[i] = - coffs[i];
        return new Polynome(newcofs);
    }
 
    /**
     * Substracts a polynome to the current polynome.
     * @param polynome the polynome to substract.
     * @return the result of the substraction.
     **/
    public Polynome substract(Polynome polynome)
    {
        int newDeg = Math.max(deg,polynome.getDeg());
 
        //sets the degree of the 2 polynomes to add
	Polynome thisTmp = new Polynome(newDeg);        
	for(int i = 0 ; i <= deg ; i++)
	    thisTmp.setCof(i,getCof(i));
 
	Polynome polTmp = new Polynome(newDeg);
	for(int i = 0 ; i <= polynome.getDeg() ; i++)
	    polTmp.setCof(i,polynome.getCof(i));
 
        //adds the polynomes
	for(int i = 0 ; i <= newDeg ; i++)
	    thisTmp.setCof(i,thisTmp.getCof(i) - polTmp.getCof(i));
        thisTmp.trim();
        return thisTmp;
    }        
 
    /**
     * Multiply the current polynome by a scalar.
     * @param scalar the scalar wich multiplies the polynome.
     * @return the resutl of the multiplication.
     **/
    public Polynome multByScal(double scalar)
    {    
        if(scalar == 0)
            return new Polynome(0);
        else if(scalar == 1)
            return clone();
        else            
        {
            Polynome thisTmp = clone();
            for(int i = 0 ; i <= deg ; i ++)
                thisTmp.setCof(i,getCof(i) * scalar);
            return thisTmp;
        }
    }        
 
    /**
     * Multiplies a polynome to the current polynome
     * @param polynome the polynome to multiply.
     * @return the result of the multiplication.
     **/
    public Polynome multiply(Polynome polynome)
    {
        Polynome pol = polynome.clone();
	int newDeg = deg + pol.getDeg();
 
	Polynome thisTmp = new Polynome(newDeg);        
	for(int i = 0 ; i <= deg ; i++)
	    thisTmp.setCof(i, getCof(i));        
 
	Polynome polTmp = new Polynome(newDeg);
	for(int i = 0 ; i <= pol.getDeg() ; i++)
	    polTmp.setCof(i,pol.getCof(i));
 
        int count = newDeg;
        while(count >= 0)
        {
            double sum = 0;
            for(int i = 0 ; i <= newDeg ; i++)
                for(int j = 0 ; j <= newDeg ; j++)
                    if(i + j == count)
                        sum += thisTmp.getCof(i) * polTmp.getCof(j);
            thisTmp.setCof(count,sum);
            count--;            
        }
 
        return thisTmp;
    }        
 
}

La classe Vector sera fournie dans les post suivants.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 
/*
 * RNFunction.java
 *
 * Created on 28 novembre 2006, 10:23
 *
 */
 
package org.umh.math.analysis;
 
import org.umh.math.linear.Vector;
 
/**
 * A function is represented by its evaluation.
 * @author Absil Romain
 */
public interface RNFunction
{
    /**
     * Evaluate the function on the given argument.
     * @return the evaluation of the function on the given argument.
     * @param v the vector you want to function to be evaluated on.
     */
    public Vector eval(Vector v);
}

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 
/*
 * Sequence.java
 *
 * Created on 30 novembre 2006, 9:38
 *
 */
 
package org.umh.math.analysis;
 
/**
 * A seuquence if a function from N to R with is represented by its evaluation.
 * @author Absil Romain
 */
public interface Sequence
{
    /**
     * Returns the sequence evaluated on the natural n.
     * @param n the natural you want the sequence to be evaluated on.
     * @return  the sequence evaluated on the natural n.
     */
    public double eval(int n);
}

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
 
/*
 * Sequence2D.java
 *
 * Created on 30 novembre 2006, 9:38
 *
 */
 
package org.umh.math.analysis;
 
import java.awt.geom.Point2D;
 
/**
 * This clas models sequences from N to R² wich are rerpesented by 
 * their evaluations.
 * @author Absil Romain
 */
public class Sequence2D
{
    private Sequence fx;
    private Sequence fy;
 
    /**
     * Construcs a new Sequence2D with the sequences defining its evaluation.
     * @param fx the sequence giving the x coordinate of the evaluations.
     * @param fy the sequence giving the y coordinate of the evaluations.
     **/
    public Sequence2D(Sequence fx, Sequence fy)
    {
        this.setFx(fx);
        this.setFy(fy);
    }
 
    /**
     * Evaluates the sequence to the given parameter.
     * @param t the parameter you want the sequence to be evaluated on.
     * @return the evaluation of the sequence to the given parameter.
     **/
    public Point2D.Double eval(int t)
    {
        return new Point2D.Double(getFx().eval(t), getFy().eval(t));
    }
 
    /**
     * Returns the sequence giving the x coordinate of the evaluations.
     * @return the sequence giving the x coordinate of the evaluations.
     **/
    public Sequence getFx()
    {
        return fx;
    }
 
    /**
     * Sets the sequence giving the x coordinate of the evaluations.
     * @param fx the sequence to set.
     **/
    public void setFx(Sequence fx)
    {
        this.fx = fx;
    }
 
    /**
     * Returns the sequence giving the y coordinate of the evaluations.
     * @return the sequence the y coordinate of the evaluations.
     **/
    public Sequence getFy()
    {
        return fy;
    }
 
    /**
     * Sets the sequence giving the y coordinate of the evaluations.
     * @param fy the sequence to set.
     **/
    public void setFy(Sequence fy)
    {
        this.fy = fy;
    }
 
 
}

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 
/*
 * Sequence_1_2.java
 *
 * Created on 22 mars 2007, 15:56
 *
 */
 
package org.umh.math.analysis;
 
import java.awt.geom.Point2D;
 
/**
 * This class models sequences from N² to R.
 * @author Absil Romain
 */
public interface Sequence_2_1
{
    /**
     * Evaluates the sequence to the given parameters n and m.
     * @param n the first number you want the function to be evaluated on.
     * @param m the second number you want the function to be evaluated on.
     * @return the sequence evaluated on the given point (n, m).
     */
    public double eval(int n, int m);
}

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
 
/*
 * ParametricFunction.java
 *
 * Created on 22 décembre 2006, 13:42
 *
 */
 
package org.umh.math.analysis;
 
import org.umh.math.linear.Vector;
 
/**
 * This class models vectorial functions wich are functions from R to Rn.
 * @author Absil Romain
 */
public class VectorialFunction
{
    private Function[] functions;
 
    /**
     * Construc a new vectorial, wich form is 
     * f : R -> Rn : t -> f(t) = ( f1(t) , ... , fn(t) ).
     * The functions f1, ... , fn are function from R to R wich give the 
     * function evaluated in t.
     * <p>e.g. to construct the function f(t) = (t, 1/t, t³) call
     * @param functions the functions f1, ... , fn wich give the function 
     * evaluated in t.
     **/
    public VectorialFunction(Function ... functions)
    {
        this.functions = functions;
    }
 
    /**
     * Returns the function evaluated in the real t.
     * @return the function evaluated in the real t.
     * @param t the number you want the function to be evaluated on.
     */
    public Vector eval(double t)
    {
        Vector v = new Vector(functions.length);
        for (int i = 0; i < functions.length; i++)
            v.setComp(i,functions[i].eval(t));
        return v;
    }
 
}

Voilà pour les classes d'ordre général ne faisant partie d'aucun des domaines ci-dessus mais utilisés un peu partout

Recherche de racine

Les classes fournies ci-dessous permettent la résolution de problèmes de recherche de racine de fonctions.

Superclasse AbstractRootFinder dont tous les chercheurs de racine héritent.
Notez qu'une condition d'arrêt est nécéssaire pour la recherche de racine. Les deux conditions fournies par la librairie sont un nbr d'itérations maximum ou une précision.
Classe :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
 
/*
 * RootFinder.java
 *
 * Created on 5 octobre 2006, 8:31
 * 
 */
 
package org.umh.math.analysis.rootFinding;
 
import java.util.Iterator;
 
 
/**
 * This class models an abstract root finder that will be able to find roots of
 * function by overriding the abstracts methods in the subclasses according to
 * a root finding algorithm.
 * @see java.util.Iterator
 * @author Absil Romain
 */
public abstract class AbstractRootFinder implements Iterator<Double>
{
    /**
     * Returns the current estimation of the function's root.
     * @return the current estimation of the function's root.
     **/
    public abstract double current();
 
    /**
     * Returns the next estimation of the function's root and updates the 
     * current estimation of the root to this value.
     * @return the next estimation of the function's root and updates the 
     * current estimation of the root to this value.
     **/
    public abstract Double next();
 
    /**
     * Returns the next estimation of the function's root without updating 
     * the current estimation of the root to this value.
     * @return the next estimation of the function's root without updating 
     * the current estimation of the root to this value.
     **/
    public abstract double getNext();
 
    /**
     * Returns true if the current estimation of the root is the exact value of
     * the function's root, return false otherwise.
     * @return true if the current estimation of the root is the exact value of
     * the function's root, return false otherwise.
     **/
    public abstract boolean hasNext();
 
    /**
     * Returns a root of the function, breaking the root finding algorithm 
     * according to the given condition.
     * @return a root of the function, breaking the root finding algorithm 
     * according to the given condition.
     * @param condition the condition wich will stop the root's evaluations.
     * @see org.umh.math.analysis.rootFinding.Cauchy
     * @see org.umh.math.analysis.rootFinding.MaxIteration
     */
    public final double findRoot(BreakCondition condition)
    {
        while(!condition.breakLoop() && hasNext())
            next();
        return current();
    }
 
    /**
     * Useless methode. Just to override the method <code>remove()</code>
     * in {@link java.util.Iterator}.
     **/
    public void remove()
    {
    }
}

Conditions d'arrêt : interface général.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 
/*
 * BreakCondition.java
 *
 * Created on 3 octobre 2006, 19:12
 * 
 */
 
package org.umh.math.analysis.rootFinding;
 
/**
 * This interface models the break condition of a root finding algorithm.
 * @see org.umh.math.analysis.rootFinding.AbstractRootFinder
 * @author Absil Romain
 */
public interface BreakCondition
{
    /**
     * Returns true if the root finding algorithm must be stop, false otherwise
     * @return true if the root finding algorithm must be stop, false otherwise
     **/
    public boolean breakLoop();
}

Condition d'arrêt par précision :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
 
/*
 * Cauchy.java
 *
 * Created on 3 octobre 2006, 19:58
 *
 */
 
package org.umh.math.analysis.rootFinding;
 
/**
 * This class models the Cauchy convergence criterion, 
 * i.e. test if |xn - xm| < Epsilon.<p>
 * In the implementing root finding classes, xn = the current estimation of the
 * root and xm = the next estimation of the root. 
 * @author Absil Romain
 */
public class Cauchy implements BreakCondition
{
 
    private AbstractRootFinder finder;
    private double epsilon;
 
    /**
     * Creates a new instance of Cauchy that will be able to break root finding
     * algorithms by Cauchy criterion.
     * 
     * @see org.umh.math.analysis.rootFinding.BreakCondition
     * @see org.umh.math.analysis.rootFinding.AbstractRootFinder#findRoot(
     * BreakCondition)
     * @param finder the finder you want to break the algorithm.
     * @param epsilon the precision of the estimation.
     */
    public Cauchy(AbstractRootFinder finder, double epsilon)
    {
        this.finder = finder;
	if(epsilon <= 0)
	    throw new IllegalArgumentException("Epsilon must be positive");
        this.epsilon  = epsilon;
    }
 
    /**
     * Returns true if the root finder algorithm must be broken, 
     * return false otherwise.
     * @return true if the root finder algorithm must be broken, 
     * return false otherwise.
     **/
    public boolean breakLoop()
    {
        return Math.abs(finder.current() - finder.getNext()) < epsilon;
    }
 
    /**
     * A good precision for root finding algorithms.
     */
    public static final double EPSILON = 10E-16;
}

condition d'arrêt par itérations maximum :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
 
/*
 * MaxIteration.java
 *
 * Created on 3 octobre 2006, 20:05
 * 
 */
 
package org.umh.math.analysis.rootFinding;
 
/**
 * This class models a break condition of root finding algorithm by setting 
 * a maximum number of iterations of the algorithm.
 * @author Absil Romain
 */
public class MaxIteration implements BreakCondition
{
 
    private int n;
    private int nbr;
 
    /**
     * 
     * Creates a new instance of MaxIteration that will be able to break root 
     * finding algorithms by setting a maximum number of iterations of the 
     * algorithm.
     * 
     * @see org.umh.math.analysis.rootFinding.BreakCondition
     * @see org.umh.math.analysis.rootFinding.AbstractRootFinder#findRoot(
     * BreakCondition )
     * @param nbr the maximum number of iterations of the root finding 
     * algorithm.
     */
    public MaxIteration(int nbr)
    {
        this.n = 0;
	if(nbr <= 0)
	    throw new IllegalArgumentException("The number of iterations " + 
					       "must be positive");
        this.nbr = nbr;
    }
 
    /**
     * Returns true if the root finder algorithm must be broken, 
     * return false otherwise.
     * @return true if the root finder algorithm must be broken, 
     * return false otherwise.
     **/
    public boolean breakLoop()
    {
        n++;
        return n > nbr;
    }
 
}

Algorithmes de recherche de racine

Voci ci-dessous 5 algorithmes de recherche de racine étandant la classe AbstractRootFinder :

Bissection :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
 
/*
 * Bissection.java
 *
 * Created on 3 octobre 2006, 19:09
 * 
 */
 
package org.umh.math.analysis.rootFinding;
 
import org.umh.math.analysis.Function;
 
 
/**
 * This class models a Bissection object that will be able to find roots of 
 * functions by the bissection algorithm.
 * @author Absil Romain
 * @see org.umh.math.analysis.Function
 * @see org.umh.math.analysis.rootFinding.AbstractRootFinder
 */
public class Bissection extends AbstractRootFinder
{
 
    private Function fct;
 
    private double b;
    private double a;
    private double x;
    private boolean hasNext;
    private double fx;
 
    /**
     * 
     * Creates a new instance of Bissection, wich will be able to find a root 
     * of the function in the interval ]a , b[ by the algorithm of bissection.
     * <br>(Precondition : f(a) * f(b) < 0)
     * @see org.umh.math.analysis.Function#eval(double a)
     * @param fct the function you want to find roots.
     * @param a a born of the interval.
     * @param b the other born of the interval.
     */
    public Bissection(Function fct, double a, double b)
    {   
        double fa = fct.eval(a);
        double fb = fct.eval(b);
        if(fa * fb > 0) //teste l'hypothèse de l'algorithme : f(a) * f(b) < 0
            throw new IllegalArgumentException("fct.eval(a) * fct.eval(b) " + 
					       "must be negative");
        if(fct.eval(a) > 0) //teste si f(a) < 0 et f(b) > 0, sinon switch(a,b)
        {
            double tmp = a;
            a = b;
            b = tmp;
        }
 
        this.fct = fct; //fonction à étudier
        this.a = a; //borne de l'intervalle
        this.b = b; //borne de l'intervalle
        this.x = (a + b) / 2; //estimation de la racine
	this.fx = fct.eval(x);
        hasNext = fx == 0 ? false : true;//on doit continuer à estimer?
    }
 
    /**
     * Returns the current estimation of the function's root.
     * @return the current estimation of the function's root.
     **/
    public double current()
    {
        return x;
    }
 
    /**
     * Returns the next estimation of the function's root and updates the 
     * current estimation of the root to this value.
     * @return the next estimation of the function's root and updates the 
     * current estimation of the root to this value.
     **/
    public Double next()
    {
        if(fx > 0)
            b = x;
        else if(fx < 0)
            a = x;
        else
        {
            hasNext = false;
            return x;
        }	
        x = (a + b) / 2;
	fx = fct.eval(x);
	return x;
    }
 
    /**
     * Returns the next estimation of the function's root without updating 
     * the current estimation of the root to this value.
     * @return the next estimation of the function's root without updating 
     * the current estimation of the root to this value.
     **/
    public double getNext()
    {
        double bTmp = b;
        double aTmp = a;
        if(fx > 0)
            bTmp = x;
        else if(fx < 0)
            aTmp = x;
        else
            return x;
        return x = (aTmp + bTmp) / 2;
    }
 
    /**
     * Returns true if the current estimation of the root is the exact value of
     * the function's root in the interval ]a , b[, return false otherwise.
     * @return true if the current estimation of the root is the exact value of
     * the function's root in the interval ]a , b[, return false otherwise.
     **/
    public boolean hasNext()
    {
        return hasNext;
    }       
 
}

Fausse position :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
 
/*
 * FalsePotision.java
 *
 * Created on 3 octobre 2006, 19:09
 * 
 */
 
package org.umh.math.analysis.rootFinding;
 
import org.umh.math.analysis.Function;
 
 
/**
 * This class models a FalsePosition object that will be able to find roots of 
 * functions by the false position algorithm.
 * @author Absil Romain
 * @see org.umh.math.analysis.Function
 * @see org.umh.math.analysis.rootFinding.AbstractRootFinder
 */
public class FalsePosition extends AbstractRootFinder
{
 
    private Function fct;
 
    private double b;
    private double a;
    private double x;
    private boolean hasNext;
    private double fx;
 
    /**
     * 
     * Creates a new instance of FalsePosition, wich will be able to find a 
     * root
     * of the function fct in the interval ]a , b[ by the algorithm of false 
     * position.
     * <br>(Precondition : f(a) * f(b) < 0)
     * 
     * @see org.umh.math.analysis.Function#eval(double a)
     * @param fct the function you want to find roots.
     * @param a a born of the interval.
     * @param b the other born of the interval.
     */
    public FalsePosition(Function fct, double a, double b)
    {
        double fa = fct.eval(a);
        double fb = fct.eval(b);
        if(fa * fb > 0) //teste l'hypothèse de l'algorithme : f(a) * f(b) < 0
            throw new IllegalArgumentException("fct.eval(a) * fct.eval(b) < " +
                    "0 must return true");
        if(fa > 0) //teste si f(a) < 0 et f(b) > 0, sinon switch(a,b)
        {
            double tmp = a;
            a = b;
            b = tmp;
        }
 
        this.fct = fct; //fonction à étudier
        this.a = a; //borne de l'intervalle
        this.b = b; //borne de l'intervalle
        this.x = a - ( (b - a) / (fb - fa) ) * fa; //estimation de la racine
	this.fx = fct.eval(x);
        hasNext = fx == 0 ? false : true;//on doit continuer à estimer?
    }
 
    /**
     * Returns the current estimation of the function's root.
     * @return the current estimation of the function's root.
     **/
    public double current()
    {
        return x;
    }
 
    /**
     * Returns the next estimation of the function's root and updates the 
     * current estimation of the root to this value.
     * @return the next estimation of the function's root and updates the 
     * current estimation of the root to this value.
     **/
    public Double next()
    {
        if(fx > 0)
            b = x;
        else if(fx < 0)
            a = x;
        else
        {
            hasNext = false;
            return x;
        }
        double fa = fct.eval(a);
        x = a - ( (b - a) / (fct.eval(b) - fa) ) * fa;
	fx = fct.eval(x);
	return x;
    }
 
    /**
     * Returns the next estimation of the function's root without updating 
     * the current estimation of the root to this value.
     * @return the next estimation of the function's root without updating 
     * the current estimation of the root to this value.
     **/
    public double getNext()
    {
        double bTmp = b;
        double aTmp = a;
        if(fx > 0)
            bTmp = x;
        else if(fx < 0)
            aTmp = x;
        else
            return x;
        double fa = fct.eval(a);
        return x = a - ( (b - a) / (fct.eval(a) - fa) ) * fa;
    }
 
    /**
     * Returns true if the current estimation of the root is the exact value of
     * the function's root in the interval ]a , b[, return false otherwise.
     * @return true if the current estimation of the root is the exact value of
     * the function's root in the interval ]a , b[, return false otherwise.
     **/
    public boolean hasNext()
    {
        return hasNext;
    }
 
}

Sécante :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
 
/*
 * Secant.java
 *
 * Created on 11 octobre 2006, 09:41
 * 
 */
 
package org.umh.math.analysis.rootFinding;
 
import org.umh.math.analysis.Function;
 
 
/**
 * This class models a Secant object that will be able to find roots of 
 * functions by the secant algorithm.
 * @author Absil Romain
 * @see org.umh.math.analysis.Function
 * @see org.umh.math.analysis.rootFinding.AbstractRootFinder
 */
public class Secant extends AbstractRootFinder
{
 
    private Function fct;
    private double x;
    private double lastX;
    private boolean hasNext;
 
    /** 
     * Creates a new instance of Secant, wich will be able to find a root of
     * the function fct in the interval ]a , b[ by the algorithm of secant.
     * @param fct the function you want to find roots.
     * @param a a born of the interval.
     * @param b the other born of the interval.
     * @see org.umh.math.analysis.Function#eval(double a)
     */
    public Secant(Function fct, double a, double b)
    {
        this.fct = fct;
        this.lastX = a;
        this.x = b;
        hasNext = fct.eval(x) == 0 ? false : true;//on doit continuer à estimer
    }
 
    /**
     * Returns the current estimation of the function's root.
     * @return the current estimation of the function's root.
     **/
    public double current()
    {
        return x;
    }
 
    /**
     * Returns the next estimation of the function's root and updates the 
     * current estimation of the root to this value.
     * @return the next estimation of the function's root and updates the 
     * current estimation of the root to this value.
     **/
    public Double next()
    {
        double fx = fct.eval(x);
        double newx = x - ( (x - lastX) / (fx - fct.eval(lastX)) ) * fx;
        lastX = x;
        x = newx;
        hasNext = fct.eval(x) == 0 ? false : true;
        return x;
    }
 
    /**
     * Returns the next estimation of the function's root without updating 
     * the current estimation of the root to this value.
     * @return the next estimation of the function's root without updating 
     * the current estimation of the root to this value.
     **/
    public double getNext()
    {
        double fx = fct.eval(x);
        return x - ( (x - lastX) / (fx - fct.eval(lastX)) ) * fx;
    }
 
    /**
     * Returns true if the current estimation of the root is the exact value of
     * the function's root in the interval ]a , b[, return false otherwise.
     * @return true if the current estimation of the root is the exact value of
     * the function's root in the interval ]a , b[, return false otherwise.
     **/
    public boolean hasNext()
    {
        return hasNext;
    }
 
}

Newton :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
 
/*
 * Newton.java
 *
 * Created on 13 octobre 2006, 08:55
 * 
 */
 
package org.umh.math.analysis.rootFinding;
 
import org.umh.math.analysis.Function;
 
 
/**
 * This class models a Newton object that will be able to find roots of 
 * functions by the Newton algorithm.
 * @see org.umh.math.analysis.Function
 * @see org.umh.math.analysis.rootFinding.AbstractRootFinder
 * @author Absil Romain
 */
public class Newton extends AbstractRootFinder
{
 
    private Function fct;
    private Function deriv;
    private double x;
    private boolean hasNext;
    private double fx;
 
    /** 
     * Creates a new instance of Newton, wich will be able to find a root of
     * the function fct near the point a by the algorithm of Newton.<br>
     * (Precondition : the derivative of the function must not be null beetween
     * the point and the root, otherwise the method 
     * {@link org.umh.math.analysis.rootFinding.AbstractRootFinder#findRoot(
     * BreakCondition )} is expected to return {@link java.lang.Double#NaN}) 
     * @param fct the function you want to find roots.
     * @param deriv the derivative of the function fct.
     * @param a the point near the root.
     */
    public Newton(Function fct, Function deriv, double a)
    {
        this.fct = fct;
        this.deriv = deriv;
        this.x = a;
	this.fx = fct.eval(x);
        hasNext = fx == 0 ? false : true;//on doit continuer à estimer?
    }
 
    /**
     * Returns the current estimation of the function's root.
     * @return the current estimation of the function's root.
     **/
    public double current()
    {
        return x;
    }
 
    /**
     * Returns the next estimation of the function's root and updates the 
     * current estimation of the root to this value.
     * @return the next estimation of the function's root and updates the 
     * current estimation of the root to this value.
     **/
    public Double next()
    {
        x = x - fx / deriv.eval(x);
	fx = fct.eval(x);
        if(fx == 0)
            hasNext = false;
        return x;
    }
 
    /**
     * Returns the next estimation of the function's root without updating 
     * the current estimation of the root to this value.
     * @return the next estimation of the function's root without updating 
     * the current estimation of the root to this value.
     **/
    public double getNext()
    {
        return x - fct.eval(x) / deriv.eval(x);
    }
 
    /**
     * Returns true if the current estimation of the root is the exact value of
     * the function's root near the point a, return false otherwise.
     * @return true if the current estimation of the root is the exact value of
     * the function's root near the point a, return false otherwise.
     **/
    public boolean hasNext()
    {
        return hasNext;
    }
 
}

Point fixe :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
 
/*
 * FixPoint.java
 *
 * Created on 13 octobre 2006, 08:45
 *
 */
 
package org.umh.math.analysis.rootFinding;
 
import org.umh.math.analysis.Function;
 
 
/**
 * Thi class models a FixPoint object that will be able a roots of a function 
 * on the line y = x by the algorithm of fix point.<br>
 * (Precondition : The absolute value of the derivative calculated in the fix 
 * point is expected to be < 1, otherwise the method may not converge)
 * <br>Warning : do not use {@link org.umh.math.analysis.rootFinding.Cauchy} for 
 * root finding or the algorithm may run in an infinite loop in case of 
 * repulsive fix point, 
 * i.e.{@link org.umh.math.analysis.rootFinding.AbstractRootFinder#findRoot(
 * BreakCondition)} may not converge to the root.
 * @author Absil Romain
 * @see org.umh.math.analysis.Function
 * @see org.umh.math.analysis.rootFinding.AbstractRootFinder
 */
public class FixPoint extends AbstractRootFinder
{
 
    private Function fct;
    private double x;
    private double next;
    private boolean hasNext;
 
    /** 
     * Creates a new instance of FixPoint, wich will be able to find a root of
     * the function fct on the line y = x near the point a by the algorithm of 
     * fix point.
     * <br>Warning : do not use {@link org.umh.math.analysis.rootFinding.Cauchy} 
     * for root finding or the algorithm may run in an infinite loop in case of 
     * repulsive fix point.
     * @param fct the function you want to find roots.
     * @param a the point near the root.
     */
    public FixPoint(Function fct, double a)
    {        
        this.fct = fct;
        this.x = a;        
        hasNext = fct.eval(x) == 0 ? false : true;//on doit continuer à estimer
    }
 
    /**
     * Returns the current estimation of the function's root.
     * @return the current estimation of the function's root.
     **/
    public double current()
    {
        return x;
    }
 
    /**
     * Returns the next estimation of the function's root and updates the 
     * current estimation of the root to this value.
     * @return the next estimation of the function's root and updates the 
     * current estimation of the root to this value.
     **/
    public Double next()
    {
        x = fct.eval(x);
        if(fct.eval(x) == 0)
            hasNext = false;
        return x;
    }
 
    /**
     * Returns the next estimation of the function's root without updating 
     * the current estimation of the root to this value.
     * @return the next estimation of the function's root without updating 
     * the current estimation of the root to this value.
     **/
    public double getNext()
    {
        return fct.eval(x);
    }
 
    /**
     * Returns true if the current estimation of the root is the exact value of
     * the function's root near the point a, return false otherwise.
     * @return true if the current estimation of the root is the exact value of
     * the function's root near the point a, return false otherwise.
     **/
    public boolean hasNext()
    {
        return hasNext;
    }
 
}

Résolution de systèmes linéaires

Voici maintenant plusieurs classes utiles à la résolution de systèmes linéaires et en algèbre linéaire. Je les poste en plusisurs fois pq elles sont longue et ça fait planter les posts...

Classe Vector :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
 
/*
 * Vector.java
 *
 * Created on 22 mars 2007, 15:59
 *
 */
 
package org.umh.math.linear;
 
import java.awt.geom.Point2D;
 
/**
 * This class models a Vector object that will be able to
 * @author Absil Romain
 */
public class Vector
{
    /**
     * Creates a null vector with n composants
     * @param dim the vector's number composants.
     **/
    public Vector(int dim)
    {
        this.dim = dim;
        comp = new double[dim];
    }
 
    /**
     * Creates a vector giving explicitly it's composants.
     * @param cofs the composants of the vector to create.
     **/
    public Vector(double ... cofs)
    {
        dim = cofs.length;
        comp = new double[cofs.length];
        for(int i = 0; i < cofs.length ; i++)
            comp[i] = cofs[i];
    }
 
    /**
     * Returns the dimension of the current vector.
     * @return the dimension of the current vector.
     **/
    public int size()
    {
        return dim;
    }
 
    /**
     * Returns the i_th composant of the current vector.
     * @return the i_th composant of the current vector.
     * @param i the index of the composant you want to return.
     */
    public double getComp(int i)
    {
        return comp[i];
    }
 
    /**
     * Returns all the composants of the current vector.
     * @return all the composants of the current vector.
     **/
    public double[] getComposants()
    {
        return comp;
    }
 
    /**
     * Sets the i_th composant of the current vector.
     * @param val the value you want to set.
     * @param i the i_th composant you want to set.
     * éparam val the value you want to put.
     */
    public void setComp(int i , double val)
    {
        comp[i] = val;
    }
 
    /**
     * Returns true if the current matrix is equals to the given object, returns
     * false otherwise.
     * @param other the object you want to know if it is equals to the current matrix.
     * @return true if the current matrix is equals to the given object, false otherwise.
     */
    public boolean equals(Object other)
    {
        if(!(other instanceof Vector))
            throw new IllegalArgumentException();
        else
        {
            Vector vect = (Vector)other;
            double[] vectComp = vect.getComposants();
            for (int i = 0; i < dim; i++)
                if(comp[i] != vectComp[i])
                    return false;
            return true;
        }
    }
 
    /**
     * Returns a copy of the current vector.
     * @return a copy of the current vector.
     **/
    public Vector clone()
    {
        return new Vector(comp);
    }
 
    /**
     * Returns the string representation of the current vector.
     * @return the string representation of the current vector.
     **/
    public String toString()
    {
        String vector = "( ";
        for (int i = 0; i < dim; i++)
        {
            if(i == dim - 1)
                vector += comp[i] + " )";
            else
                vector += comp[i] + " , ";
        }
        return vector;
    }
 
    /**
     * Parses the current vector to a point with coordinates (x, y), if the
     * current vector has exactly two composants.
     * @return the current vector parsed to a point with coordinates (x , y).
     */
    public Point2D.Double toPoint2D()
    {
        if(dim != 2)
            throw new IllegalArgumentException("The vector must have " +
                    "exactly two composants");
        return new Point2D.Double(comp[0], comp[1]);
    }
 
    /**
     * Returns the norm of the current polynome.
     * @return the norm of the current polynome.
     **/
    public double getNorm()
    {
        double sum = 0;
        for(int i = 0 ; i < dim ; i++)
            sum += comp[i] * comp[i];
        return Math.sqrt(sum);
    }
 
    /**
     * Adds a vector to the current vector.
     * <p>(Precondition : dim == vect.getDim() )
     * @param vect the vector to add.
     * @return the result of the addition.
     **/
    public Vector add(Vector vect)
    {
        if(dim != vect.size())
            throw new IllegalArgumentException();
        Vector added = vect.clone();
        for(int i = 0 ; i < dim ; i++)
            added.setComp(i, added.getComp(i) + comp[i]);
        return added;
    }
 
    /**
     * Returns the opposite of the current vector.
     * @return the opposite of the current vector.
     **/
    public Vector getOpposite()
    {
        double[] cofs = new double[comp.length];
        for (int i = 0; i < cofs.length; i++)
            cofs[i] = -comp[i];
        return new Vector(cofs);
    }
 
    /**
     * Mutliplies the curent vector by the given scalar.
     * @param scalar the scalar wich multiplies the current vector
     * @return the result of the addition.
     **/
    public Vector multByScal(double scalar)
    {
        if(scalar == 0)
            return new Vector(dim);
        else if(scalar == 1)
            return clone();
        else
        {
            Vector vector = clone();
            for(int i = 0 ; i < dim ; i++)
                vector.setComp(i,vector.getComp(i) * scalar);
            return vector;
        }
    }
 
    /**
     * Substract the given vector to the current vector.
     * <p>(Precondition : the dimensions of the two vectors must be the same)
     * @param vect the vector to substract.
     * @return the result of the substraction.
     **/
    public Vector substract(Vector vect)
    {
        if(dim != vect.size())
            throw new IllegalArgumentException();
        Vector added = vect.clone();
        for(int i = 0 ; i < dim ; i++)
            added.setComp(i, -added.getComp(i) + comp[i]);
        return added;
    }
 
    /**
     * Returns the algebric scalar product
     * <p>(Precondition : dim == vect.getDim() )
     * @param vect the vector you want to multiply.
     * @return the current vector scalar the given vector.
     **/
    public double multiply(Vector vect)
    {
        if(dim != vect.size())
            throw new IllegalArgumentException();
        double sum = 0;
        for(int i = 0 ; i < dim ; i++)
            sum += comp[i] * vect.getComp(i);
        return sum;
    }
 
    private int dim;
    private double[] comp;
}

Classe Matrice :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
 
/*
 * Matrix.java
 *
 * Created on 22 mars 2007, 16:01
 *
 */
 
package org.umh.math.linear;
 
/**
 * This class models a Matrix object that will be able to
 * @author Absil Romain
 */
public class Matrix
{
    /**
     * The rowspace of the matrix (i.e. the set of its rows view as a set of row
     * vectors.
     */
    protected Vector[] rowSpace;
    /**
     * The number of columns of the matrix.
     */
    protected int n;
    /**
     * The number of rows of the matrix.
     */
    protected int m;
 
    /**
     * Conctructs a new num matrix with m columns and n rows.
     * @param m the number of rows of the matrix.
     * @param n the number of columns of the matrix.
     **/
    public Matrix(int m, int n)
    {
        this.rowSpace = new Vector[m];
        for (int i = 0; i < rowSpace.length; i++)
            rowSpace[i] = new Vector(n);
        this.n = m;
        this.m = n;
    }
 
    /**
     * Constructs a new matrix with the given 2D array.
     * @param cofs the coefficients of the matrix.
     **/
    public Matrix(double[][] cofs)
    {
        this.n = cofs[0].length;
        this.m = cofs.length;
        this.rowSpace = new Vector[m];
 
        for (int i = 0; i < m; i++)
            rowSpace[i] = new Vector(cofs[i]);
    }
 
    /**
     * Construcs a new matrix with the rows are the given vectors.
     * @param vectors the rows of the matrix.
     **/
    public Matrix(Vector ... vectors)
    {
        this.m = vectors.length;
        this.n = vectors[0].size();
        rowSpace = new Vector[m];
        rowSpace[0] = vectors[0];
        for (int i = 1; i < m; i++)
        {
            if(vectors[i].size() != n)
                throw new IllegalArgumentException(
                        "All vectors must be the same dimension");
            rowSpace[i] = vectors[i];
        }
    }
 
    /**
     * Returns the line space of the matrix.
     * @return the line space of the matrix.
     **/
    public Vector[] getRowSpace()
    {
        return rowSpace;
    }
 
    /**
     * Returns the number of rows of the current matrix.
     * @return the number of rows of the current matrix.
     **/
    public int getRowsNumber()
    {
        return m;
    }
 
    /**
     * Returns the number of columns of the current matrix.
     * @return the number of columns of the current matrix.
     **/
    public int getColumnsNumber()
    {
        return n;
    }
 
    /**
     * Returns the element (i,j) of the current matrix.
     * @param i the row of the element to return.
     * @param j the column of the element to return.
     * @return the element (i,j) of the current matrix.
     **/
    public double get(int i, int j)
    {
        return rowSpace[i].getComp(j);
    }
 
    /**
     * Returns the row indexed by the given number of the current matrix.
     * @return the row indexed by the given number of the current matrix.
     * @param i the index of the row you want to be returned.
     */
    public Vector getRow(int i)
    {
        return rowSpace[i];
    }
 
    /**
     * Returns the column indexed by the given number of the current matrix.
     * @return the column indexed by the given number of the current matrix.
     * @param i the index of the columns you want to be returned.
     */
    public Vector getColumn(int i)
    {
        Vector vector = new Vector(m);
        for (int j = 0; j < m; j++)
            vector.setComp(j,rowSpace[j].getComp(i));
        return vector;
    }
 
    /**
     * Sets the element (i,j) to the given value.
     * @param i the row of the element to set.
     * @param j the column of the element to set.
     * @param value the element to set.
     **/
    public void set(int i, int j, double value)
    {
        rowSpace[i].setComp(j,value);
    }
 
    /**
     * Sets the line of the given matrix to the given vector.
     * @param i the line to set.
     * @param vector the nex value of the line.
     **/
    public void setRow(int i, Vector vector)
    {
        rowSpace[i] = vector;
    }
 
    /**
     * Sets the column of the given matrix to the given vector.
     * @param i the column to set.
     * @param vector the nex value of the column.
     **/
    public void setColumn(int i, Vector vector)
    {
        for (int j = 0; j < m; j++)
            rowSpace[j].setComp(i,vector.getComp(j));
    }
 
    /**
     * Extracts a submatrix from the current matrix at the elements (l,c) and the rank (height, width).
     * @param l the index of the row you want the matrix to be extracted.
     * @param c the index of the column you want the matrix to be extracted
     * @param height the height of the extracted matrix
     * @param width the width of the extracted matrix
     * @return a submatrix from the current matrix at the elements (l,c) and the rank (height, width).
     */
    public Vector[] subMatrix(int l, int c, int height, int width)
    {
        Vector[] matrix = new Vector[height];
        for(int i = 0; i < height; i++)
        {
            Vector v = new Vector(width);
            for(int j = 0; j < width; j++)
                v.setComp(j,get(l+i,c+j));
            matrix[i] = v;
        }
        return matrix;
    }
 
    /**
     * Returns true if the given object and the current matrix are equals, 
     * returns false otherwise.
     * @param other the object you want to know if it is equals to the current matrix.
     * @return true if the given object and the current matrix are equals, 
     * returns false otherwise.
     */
    public boolean equals(Object other)
    {
        if(!(other instanceof Matrix))
            return false;
        Matrix matrix = (Matrix)other;
        for (int i = 0; i < m; i++)
            if(!rowSpace[i].equals(matrix.getRow(i)))
                return false;
        return true;
    }
 
    /**
     * Returns a copy of the current matrix.
     * @return a copy of the current matrix.
     */
    public Matrix clone()
    {
        return new Matrix(rowSpace);
    }
 
    /**
     * Returns the String representation of the current matrix.
     * @return the String representation of the current matrix.
     **/
    public String toString()
    {
        String s = "";
	for(int i = 0 ; i < m ; i++)
	    for(int j = 0 ; j < n ; j++)
		{
		    s += get(i,j) + "\t";
		    if(j == n-1)
			s += "\n";
		}
	return s;
    }
 
    /**
     * Adds the given matrix to the current matrix and return the result of 
     * the addition.
     * @param matrix the matrix to add.
     * @return the result of the addition of the current matrix to the given other matrix.
     */
    public Matrix add(Matrix matrix)
    {
        if(n != matrix.getRowsNumber() || m != matrix.getColumnsNumber())
            throw new IllegalArgumentException();
        Matrix added = clone();
        for (int i = 0; i < m; i++)
            added.setRow(i,matrix.getRow(i).add(matrix.getRow(i)));
        return added;
    }        
 
    /**
     * Returns the opposite of the current matrix for additive law.
     * @return the opposite of the current matrix for additive law.
     */
    public Matrix getOpposite()
    {
        Vector[] oppositeRows = new Vector[m];
        for (int i = 0; i < oppositeRows.length; i++)
            oppositeRows[i] = rowSpace[i].getOpposite();
        return new Matrix(oppositeRows);
    }
 
    /**
     * Substracts the given matrix to the current matrix.
     * @param matrix the matrix to substract
     * @return the result of the substraction of the current matrix and the given other matrix.
     */
    public Matrix substract(Matrix matrix)
    {
        if(n != matrix.getRowsNumber() || m != matrix.getColumnsNumber())
            throw new IllegalArgumentException();
        Matrix added = clone();
        for (int i = 0; i < m; i++)
            added.setRow(i,matrix.getRow(i).substract(matrix.getRow(i)));
        return added;
    }        
 
    /**
     * Multiplies the current matrix by the given scalar, i.e. by a real number.
     * @param scalar the number you want the current matrix to be multiplied by.
     * @return the result of the multiplication of the current matrix by the given scalar.
     */
    public Matrix multByScal(double scalar)
    {
        Matrix added = clone();
        for (int i = 0; i < m; i++)
            added.setRow(i,added.getRow(i).multByScal(scalar));
        return added;
    }
 
    /**
     * Multiplies the current matrix by another matrix.
     * @param matrix the matrix you want multiply the current matrix by.
     * @return the result of the multiplication of the current matrix by the given other matrix.
     */
    public Matrix multiply(Matrix matrix)
    {                
        Matrix mult = clone();
        for (int i = 0; i < m; i++)
            for (int j = 0; j < n; j++)
                mult.set(i,j,mult.getRow(i).multiply(matrix.getColumn(j)));
        return mult;
    }        
 
    /**
     * Sets the upper part of the matrix triangular and return the factor by
     * wich the determinant must be multiplied in case of suare matrix.
     * @return the number by wich the determinant of the matrix must be 
     * multiplied by after the triangularisation.
     */
    public double setTriangularSup()
    {
        int rowIndex = 0;
        int columnIndex = 0;
        double factor = 1; //factor to return
        while(columnIndex < n && rowIndex < m)
        {
            //if the pivot = 0, swap line
            if(get(rowIndex,columnIndex) == 0)
            {
                int index = indexofFirstDownNonNullPivot(columnIndex);
                if(index != -1) //pivot !=0 found
                {
                    swapRow(rowIndex,index);
                    factor *= -1;
                }
                else //if no pivot != 0, nothing to do
                {
                    columnIndex++;
                    factor = 0; //det = 0 
                    continue; //don't make the rest of the loop
                }
            }
 
            //make null all the elements abose the pivot
            // L_i <- L_i - k * L_{pivot != 0} and k choosen to make the 
            // element i,columnIndex null
            double pivot = get(rowIndex, columnIndex); 
            for(int i = rowIndex + 1; i < m ; i++)
                setRow(i,rowSpace[i].substract(
                        rowSpace[rowIndex].multByScal(get(i,columnIndex)  
                        / pivot)));
 
            //increment line and column index to continue the algorithm
            rowIndex++;
            columnIndex++;
        }
        return factor;
    }
 
    public double setTriangularSup(Vector independant)
    {
        int rowIndex = 0;
        int columnIndex = 0;
        double factor = 1; //factor to return
        while(columnIndex < n && rowIndex < m)
        {
            //if the pivot = 0, swap line
            if(get(rowIndex,columnIndex) == 0)
            {
                int index = indexofFirstDownNonNullPivot(columnIndex);
                if(index != -1) //pivot !=0 found
                {
                    swapRow(rowIndex,index);
                    factor *= -1;
                }
                else //if no pivot != 0, nothing to do
                {
                    columnIndex++;
                    factor = 0; //det = 0 
                    continue; //don't make the rest of the loop
                }
            }
 
            //make null all the elements abose the pivot
            // L_i <- L_i - k * L_{pivot != 0} and k choosen to make the 
            // element i,columnIndex null
            double pivot = get(rowIndex, columnIndex); 
            for(int i = rowIndex + 1; i < m ; i++)
            {
                if(independant != null)
                    independant.setComp(i, independant.getComp(i) - 
                            independant.getComp(rowIndex) * 
                            get(i,columnIndex) / pivot);
 
                setRow(i,rowSpace[i].substract(
                        rowSpace[rowIndex].multByScal(get(i,columnIndex)  
                        / pivot)));
            }
 
            //increment line and column index to continue the algorithm
            rowIndex++;
            columnIndex++;
        }
        return factor;
    }
 
    /**
     * Sets the lower part of the matrix triangular and return the factor by
     * wich the determinant must be multiplied in case of suare matrix.
     * @return the number by wich the determinant of the matrix must be multiplied by after the triangularisation.
     */
    public double setTriangularInf()
    {
        int rowIndex = n - 1;
        int columnIndex = m - 1;
        double factor = 1;
        while(rowIndex >= 0 && columnIndex >= 0)
        {
            //if the pivot = 0, swap line
            if(get(rowIndex,columnIndex) == 0)
            {
                int index = indexofFirstUpNonNullPivot(columnIndex);
                if(index != -1) //pivot !=0 found
                {
                    swapRow(rowIndex,index);
                    factor *= -1;
                }
                else //if no pivot != 0, nothing to do
                {
                    columnIndex--;
                    factor = 0; //det = 0 
                    continue; //don't make the rest of the loop
                }
            }
 
            //make null all the elements abose the pivot
            // L_i <- L_i - k * L_{pivot != 0} and k choosen to make the 
            // element i,columnIndex null
            double pivot = get(rowIndex, columnIndex); 
            for(int i = rowIndex - 1; i >= 0 ; i--)
                setRow(i,rowSpace[i].substract(
                        rowSpace[rowIndex].multByScal(get(i,columnIndex)  
                        / pivot)));
 
            //increment line and column index to continue the algorithm
            rowIndex--;
            columnIndex--;
        }
 
        return factor;
    }
 
    // swaps the line i and the line j
    private void swapRow(int i, int j)
    {
        Vector v = rowSpace[i];
        rowSpace[i] = rowSpace[j];
        rowSpace[j] = v;
    }
 
    //returns the index of the first no null pivot of the matrix at the 
    //given column, returns -1 if no non null pivots at the given column 
    private int indexofFirstDownNonNullPivot(int column)
    {
        Vector vector = getColumn(column);
        for (int i = 0; i < m; i++)
            if(vector.getComp(i) == 0)
                return i;
        return -1;
    }
 
    private int indexofFirstUpNonNullPivot(int column)
    {
        Vector vector = getColumn(column);
        for (int i = m - 1; i >= 0; i--)
            if(vector.getComp(i) == 0)
                return i;
        return -1;
    }
}

Classe Matrice carrée : (methode d'inverse non implémentée, mais compte tenu des méthodes de triangularisation fournie il ne reste plus beaucoup de travail).

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
 
/*
 * SquareMatrix.java
 *
 * Created on 22 mars 2007, 16:05
 *
 */
 
package org.umh.math.linear;
 
/**
 * This class models a SquareMatrix object that will be able to
 * @author Absil Romain
 */
public class SquareMatrix extends Matrix
{
    /**
     * Creates a new instance of SquareMatrix
     * @param n the number of rows of the square matrix.
     */
    public SquareMatrix(int n)
    {
        super(n,n);
    }
 
    /**
     * Constructs a new square matrix with the specified coefficients.
     * @param cofs the coefficients of the square matrix.
     */
    public SquareMatrix(double[][] cofs)
    {
        super(cofs);
        if(cofs.length != cofs[0].length)
            throw new IllegalArgumentException("The matrix must be a square");
    }
 
    /**
     * Constructs a new square matrix with the specified vectors as rows.
     * @param vectors the rows of the square matrix.
     */
    public SquareMatrix(Vector ... vectors)
    {
        super(vectors);
        if(getRowsNumber() != getColumnsNumber())
            throw new IllegalArgumentException("The matrix must be a square");
    }
 
    /**
     * Returns the determinant of the square matrix.
     * @return the determinant of the square matrix.
     */
    public double getDet()
    {
        if(m == 1)
            return rowSpace[0].getComp(0);
        else if(m == 2)
            return rowSpace[0].getComp(0) * rowSpace[1].getComp(1) 
                - rowSpace[1].getComp(0) * rowSpace[0].getComp(1);
        else if(m == 3)
            return rowSpace[0].getComp(0) * rowSpace[1].getComp(1) * 
                    rowSpace[2].getComp(2) + rowSpace[1].getComp(0) * 
                    rowSpace[2].getComp(1) * rowSpace[0].getComp(2) + 
                    rowSpace[0].getComp(1) * rowSpace[1].getComp(2) * 
                    rowSpace[2].getComp(0) - rowSpace[2].getComp(0) * 
                    rowSpace[1].getComp(1) * rowSpace[0].getComp(2) - 
                    rowSpace[1].getComp(0) * rowSpace[0].getComp(1) * 
                    rowSpace[2].getComp(2) - rowSpace[0].getComp(0) * 
                    rowSpace[2].getComp(1) * rowSpace[1].getComp(2);
        else
        {
            //factor by wich the determinant must be multiplied
            Matrix origin = clone();
            double factor = setTriangularSup(); 
            //product of diagonal elements
            factor *= getTrace();
 
            //sets the matrix to original form
            for (int i = 0; i < rowSpace.length; i++)
                rowSpace[i] = origin.getRow(i);
 
            return factor;
        }
    }
 
    /**
     * Returns the trace of the square matrix, i.e. returns the products of the
     * diagonal elements.
     * @return the trace of the square matrix.
     */
    public double getTrace()
    {
        double trace = 1;
        for (int i = 0; i < rowSpace.length; i++)
        {
            double e = rowSpace[i].getComp(i);
            if(e == 0)
                return 0;
            else
                trace *= e;
        }
        return trace;
    }
 
    /**
     * Returns the inverse of the square matrix for the multiplicative law.
     * @return the inverse of the square matrix for the multiplicative law.
     */
    public SquareMatrix getInverse()
    {
        return new SquareMatrix(0);
    }
 
}

Résolution de systèmes linéaires

Voici enfin les classes derésolution de systèmes linéaires (seul l'algorithme de Gauss a été fourni ici, il suffira d'étednre AbstractLinearSystemSolver pour implémenter vos propres algos :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
 
/*
 * AbstractLinearSystemSolver.java
 *
 * Created on 22 mars 2007, 16:12
 *
 */
 
package org.umh.math.linear.systemsolving;
 
import org.umh.math.linear.Vector;
 
/**
 * This class models a AbstractLinearSystemSolver object that will be able to
 * solves linear system with an implemented algorithm.
 * @author Absil Romain
 */
public abstract class AbstractLinearSystemSolver
{
    /**
     * Solves the linear system on the given vector and returns the solution 
     * as a vector.
     * @param v the vector on wich you want to solve the system.
     * @return the solution of the system solved on the given vector.
     */
    public abstract Vector solve(Vector v);
}

Algo de Gauss (pivotage) :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
 
/*
 * GaussSystemSolver.java
 *
 * Created on 22 mars 2007, 16:17
 *
 */
 
package org.umh.math.linear.systemsolving;
 
import org.umh.math.linear.Matrix;
import org.umh.math.linear.Vector;
 
/**
 * This class models a GaussSystemSolver object that will be able to solve
 * systems with the Gauss algorithm.
 * @author Absil Romain
 */
public class GaussSystemSolver extends AbstractLinearSystemSolver
{
    private Matrix matrix;
 
    /**
     * Construcs a new linear system with the given matrix.
     * @param matrix the system under matricial form.
     */
    public GaussSystemSolver(Matrix matrix)
    {
        this.matrix = matrix;
    }
 
    //résout système triangulaire sur base d'un vecteur v 
    //(résout A . x = v avec A triang. )
    /**
     * Solves a linear system that is under high triangular form on the given vector.
     * @param vector the vector on wich you want to solve the system.
     * @return the solution of the system solved on the given vector.
     */
    public Vector solveTriangularSystem(Vector vector)
    {
        if(vector.size() != matrix.getRowsNumber())
            throw new IllegalArgumentException("Invalid vector");
        Vector rep = new Vector(vector.size());
        rep.setComp(vector.size() - 1,matrix.get(
                matrix.getRowsNumber()-1,matrix.getColumnsNumber()-1));
        for(int i = vector.size() - 1 ; i >= 0 ; i--)
        {
            double comp = 1/matrix.get(i,i); //1/aii
            double val = vector.getComp(i); //bi
            double sum = 0;
            for(int j = i+1 ; j < vector.size() ; j++)
                sum += matrix.get(i,j) * rep.getComp(j);
            val -= sum;
            comp *= val;
            rep.setComp(i,comp);
        }
 
        return rep;
    }
 
    /**
     * Solves the linear system on the given vector and returns the solution 
     * as a vector.
     * @param v the vector on wich you want to solve the system.
     * @return the solution of the system solved on the given vector.
     */
    public Vector solve(Vector v)
    {
        matrix.setTriangularSup(v);
        return this.solveTriangularSystem(v);
    }
 
}

Interpolation

Voici maintenant des classes d'interpolation et de moindres carrés (seul les algorithmes polynomiaux sont fournis ici, d'autres seront peut-etre à paraître sur les exponentielles).

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
 
/*
 * Interpolator.java
 *
 * Created on 16 décembre 2006, 20:52
 *
 */
 
package org.umh.math.analysis.interpolator;
 
import org.umh.math.analysis.Function;
 
/**
 * An interpolator is represented by an interpolation method.
 * @author Absil Romain, Badibanga Christian, Minacapelli Joffrey.
 */
public interface Interpolator
{
    /**
     * Returns an interpolated function.
     * @return an interpolated function.
     **/
    public Function interpolate();
}

Alorithme d'interpolation des différences divisées :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
 
/*
 * DividedDifferences.java
 *
 * Created on 16 décembre 2006, 20:54
 *
 */
 
package org.umh.math.analysis.interpolator;
 
import java.awt.geom.Point2D;
import org.umh.math.analysis.Function;
import org.umh.math.analysis.Polynome;
 
/**
 * This class models a DividedDifferences object that will be able to
 * @author Absil Romain, Badibanga Christian, Minacapelli Joffrey.
 */
public class DividedDifferences implements Interpolator
{
    private double[] xi;
    private double[] yi;
 
    /**
     * Constructs a new DividedDifferences object that will be able to
     * approximate the given function on the specified nodes by divided
     * differences interpolation.
     * @param f the function to be interpolated
     * @param xi the points you want the function to be interpolated
     */
    public DividedDifferences(Function f, double[] xi)
    {
        this.xi = xi;
        this.yi = new double[xi.length];
        for(int i = 0; i < xi.length; i++)
            yi[i] = f.eval(xi[i]);
 
        calculateDividedDifferences();
    }
 
    /**
     * Conctructs a new DividedDifferences object that will be able to 
     * approximatethe given function with n +1 nodes on the interval [a, b].
     * @param f the function to approximate.
     * @param n the number of nodes you want - 1.
     * @param a a born of the interval.
     * @param b the other born of the interval.
     **/
    public DividedDifferences(Function f, int n, double a, double b)
    {
        this.xi = new double[n];
        this.yi = new double[n];
        double step = ( b - a) / (n - 1);
 
        xi[0] = a;
        yi[0] = f.eval(a);
        xi[n-1] = b;
        yi[n-1] = f.eval(b);
 
        double x = a;//abscisses des points
        for(int i = 1; i < n-1; i++)
        {
            x += step;
            xi[i] = x;
            yi[i] = f.eval(x);
        }
 
        calculateDividedDifferences();
    }
 
    /**
     * Conctructs a new DividedDifferences object that will be able to 
     * interpolate the given datas.
     * @param pts the datas to interpolate.
     **/
    public DividedDifferences(Point2D.Double ... pts)
    {
        this.xi = new double[pts.length];
        this.yi = new double[pts.length];
 
        for (int i = 0; i < pts.length; i++)
        {
            xi[i] = pts[i].getX();
            yi[i] = pts[i].getY();
        }
 
        calculateDividedDifferences();
    }
 
    private void calculateDividedDifferences()
    {
        for(int i = 1; i < yi.length; i++)
            for(int j = yi.length-1; j >= i; j--)
                yi[j] =  ((yi[j-1] - yi[j]) / ( xi[j-i] - xi[j]));
    }
 
    /**
     * Returns the divided differences of the current object.
     * @return the divided differences of the current object.
     **/
    public double[] getDividedDifferences()
    {
        return yi;
    }
 
    /**
     * Returns the interpolation by divided differences.
     * @return the interpolation by divided differences.
     **/
    public Function interpolate()
    {
        return new Polynome(yi);
    }
 
}

Régression linéaire :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
 
/*
 * LinearRegression.java
 *
 * Created on 16 décembre 2006, 21:02
 *
 */
 
package org.umh.math.analysis.interpolator;
 
import java.awt.geom.Point2D;
import org.umh.math.analysis.Function;
 
/**
 * This class models a LinearRegression object that will be able to interpolate
 * datas by linear regresion.
 * @author Absil Romain, Minacapelli Joffrey.
 */
public class LinearRegression implements Interpolator
{
    private double[] xi;
    private double[] yi;
 
    /**
     * Construcs a new LinearRegresion object that will be able to interpolate
     * the given datas by linear regression.
     * @param pts the datas to interpolate.
     **/
    public LinearRegression(Point2D.Double ... pts)
    {
        this.xi = new double[pts.length];
        this.yi = new double[pts.length];
 
        for (int i = 0; i < pts.length; i++)
        {
            xi[i] = pts[i].x;
            yi[i] = pts[i].y;
        }
    }
 
    /**
     * Returns the interpolation by linear regression.
     * @return the interpolation by linear regression.
     **/
    public Function interpolate()
    {
        double s_xy = 0 , s_x = 0 , s_y = 0 , s_x2 = 0;
        int n = xi.length;
 
        for(int i = 0 ; i < n ; i++)
        {
            s_xy = s_xy + xi[i] * yi[i];
            s_x = s_x + xi[i];
            s_y = s_y + yi[i];
            s_x2 = s_x2 + xi[i] * xi[i];
        }
        final double a = ((n * s_xy) - (s_x * s_y)) / 
                ((n * s_x2) - (s_x * s_x));
        final double b = ((s_x2 * s_y) - (s_x * s_xy)) / 
                ((n * s_x2) - (s_x * s_x));
 
        Function f = new Function()
        {
            public double eval(double x)
            {
                return a * x + b;
            }
        };
 
        return f;
    }
 
}

Moindres carrés polynomiaux :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
 
/*
 * PolynomialIntegrator.java
 *
 * Created on 30 janvier 2007, 15:40
 *
 */
 
package org.umh.math.analysis.interpolator;
 
import java.awt.geom.Point2D;
import org.umh.math.analysis.Function;
import org.umh.math.analysis.Polynome;
import org.umh.math.linear.Matrix;
import org.umh.math.linear.Vector;
import org.umh.math.linear.systemsolving.GaussSystemSolver;
import org.umh.math.utilities.MathUtil;
 
/**
 * This class models a PolynomialIntegrator object that will be able to
 * interpole datas by applying polynomial least squares to them.
 * @author Absil Romain
 */
public class PolynomialLeastSquares implements Interpolator
{
    private double[] xi;
    private double[] yi;
    private int degree;
 
    /**
     * Construcs a new PolynomialLeastSquares object that will be able to 
     * interpolate the given datas by polynomial least squares.
     * @param degree the degree of the interpolator polynome.
     * @param pts the datas to interpolate.
     **/
    public PolynomialLeastSquares(int degree, Point2D.Double ... pts)
    {
        this.xi = new double[pts.length];
        this.yi = new double[pts.length];
        this.degree = degree;
        for(int i = 0 ; i < pts.length ; i++)
        {
            xi[i] = pts[i].x;
            yi[i] = pts[i].y;
        }
    }
 
    /**
     * Returns the interpolation by polynomial least squares.
     * @return the interpolation by polynomial least squares.
     **/
    public Function interpolate()
    {
        //Build system A x = b (with x = a_N, ... , a_0)
 
        double xpown = 0;
        for (int i = 0; i < xi.length; i++)
            xpown += MathUtil.pow(xi[i], degree);
 
        Vector v = new Vector(degree + 1); //vector b
        Matrix matrix = new Matrix(degree + 1, degree + 1);//matrix A
        for (int i = 0; i < degree + 1; i++)
        {
            for (int j = 0; /**j < degree + 1*/j <= i; j++)
            {
                //if(i == j) //diagonal elements == sum x_i^n
                  //  matrix.set(i,j,xpown);
 
                //else
                //{
                    double x = 0;
                    double vi = 0;
                    for (int k = 0; k < xi.length; k++)
                    {
                        x += MathUtil.pow(xi[k], 2*degree - i - j);
                        vi += MathUtil.pow(xi[k], degree - i) * yi[k];
                    }
 
                    matrix.set(i,j,x); 
                    if(i != j)//symetric matrix
                        matrix.set(j,i,x); 
 
                    v.setComp(i, vi);
                //}
 
            }
 
        }    
 
        GaussSystemSolver solver = new GaussSystemSolver(matrix);
        Vector sol = solver.solve(v);
 
        return new Polynome(sol.getComposants());
    }
 
}

Et enfin voici pour le final des algorithmes utiles à l'intération, principalement l'algorithme de Runge-Kutta, à n dimensions (choisir RungeKutta2D pour le RungeKutta classique à une dimension, et RK4 pour le RungeKutta général).

Interface d'intégration à une dimension

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 
/*
 * Integrator.java
 *
 * Created on 16 décembre 2006, 20:32
 *
 */
 
package org.umh.math.analysis.integrator;
 
import java.awt.geom.Point2D;
 
/**
 * An integrator is represented by an integration method.
 * @author Absil Romain.
 */
public interface Integrator2D
{
    /**
     * Returns an integrated function on some points.
     * @param n the number of points you want the function to be integrated.
     * @return an integrated function on some points.
     */
    public Point2D.Double[] integrate(int n);
}

Runge-Kutta à une dimension :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
 
/*
 * RungeKutta.java
 *
 * Created on 16 décembre 2006, 22:53
 *
 */
 
package org.umh.math.analysis.integrator;
 
import java.awt.geom.Point2D;
import org.umh.math.analysis.Function_2_1;
 
/**
 * This class models a RungeKutta object that will be able to integrate 
 * functions from R to R under the Runge-Kuta algoritmh of order 4.
 * @author Absil Romain.
 */
public class RungeKutta2D implements Integrator2D
{
    private double u0;
    private double t1;
    private double t0;
    private Function_2_1 f;
 
    /**
     * Construcs a new RungeKutta2D object with the specified function 
     * to integrate on the interval [t0, t1] and wich value on t0 is u0.
     * <br>i.e. Considering the following Cauchy problem : 
     * <br> d u(t) = f(t, u(t)) and f(t0) = u0
     * <br> the algotihm will return the function f integrated on the 
     * interval [t0, t1].
     * @param f the function to integrate.
     * @param t0 the lower born of the interval.
     * @param t1 the higher born of the interval.
     * @param u0 the function evaluated in t0.
     **/
    public RungeKutta2D(Function_2_1 f, double t0, double t1, double u0)
    {
        this.f = f;
        this.t0 = t0;
        this.t1 = t1;
        this.u0 = u0;
    }        
 
    /**
     * Returns the function integrated on n nodes in the interval [t0, t1].
     * The return type is some points of the integrated function.
     * @param n the number of nodes on wich you want the function to be 
     * integrated.
     * @return the function integrated on n nodes in the interval [t0, t1].
     **/
    public Point2D.Double[] integrate(int n)
    {
        Point2D.Double[] pts = new Point2D.Double[n];
        double ui = u0;
        pts[0] = new Point2D.Double(t0,u0);
        double h = (t1 - t0) / (n - 1);
 
        for (int i = 1; i < pts.length; i++)
        {
            double ti = t0 + i * h;
            double k1 = h * f.eval(ti, ui);
            double k2 = h * f.eval(ti + h/2 , ui + k1/2);
            double k3 = h * f.eval(ti + h/2 , ui+ k2/2);
            double k4 = h * f.eval(ti + h, ui + k3);
 
            ui += (1/6.) * (k1 + 2*k2 + 2*k3 + k4);
            pts[i] = new Point2D.Double(ti, ui);
        }
 
        return pts;
    }
 
}

Interface d'intégration à n dimensions :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 
/*
 * RNIntegrator.java
 *
 * Created on 20 décembre 2006, 10:49
 *
 */
 
package org.umh.math.analysis.integrator;
 
/**
 * An integrator is represented by an integration method.
 * @author Absil Romain.
 */
public interface RNIntegrator
{
    /**
     * Returns an integrated function on some points.
     * @return an integrated function on some points.
     * @param n the number of points you want the functio n to be integrated.
     */
    public VectorialPoint[] integrate(int n);
}

Classe VectorialPoint sevrant à stocker les couples t_i, u_i retournés par Runge-Kutta à n dimansions :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
 
/*
 * VectorialPoint.java
 *
 * Created on 20 décembre 2006, 10:49
 *
 */
 
package org.umh.math.analysis.integrator;
 
import org.umh.math.linear.Vector;
 
/**
 * This class models a VectorialPoint object that will be able to store an 
 * real number and vector. This is quite useful for functions from R to Rn,
 * (see {@link org.umh.math.analysis.integrator.RKFunction} ) because points of the graph of 
 * these functions are the couple (x, v) where x is real and v is a vector.
 * @author Absil Romain, Badibanga Christian, Minacapelli Joffrey.
 */
public class VectorialPoint
{
 
    private double x;
    private Vector y;
 
    /**
     * Conctructs a new vectorial point (x, y) with the specified real number x
     * and vector y
     * @param x the real number.
     * @param y the vector.
     **/
    public VectorialPoint(double x, Vector y)
    {
        this.setX(x);
        this.setY(y);
    }
 
    /**
     * Returns the real number x of the couple (x, y) where is a real number 
     * and y a vector.
     * @return the real number of the vectorial point.
     **/
    public double getX()
    {
        return x;
    }
 
    /**
     * Sets the real number x of the couple (x, y) where is a real number 
     * and y a vector.
     * @param x the real number of the vectorial point to set.
     **/
    public void setX(double x)
    {
        this.x = x;
    }
 
    /**
     * Returns the vector x of the couple (x, y) where is a real number 
     * and y a vector.
     * @return the vector of the vectorial point.
     **/
    public Vector getY()
    {
        return y;
    }
 
    /**
     * Sets the real number x of the couple (x, y) where is a real number 
     * and y a vector.
     * @param y the real number of the vectorial point to set.
     **/
    public void setY(Vector y)
    {
        this.y = y;
    }
 
    /**
     * Returns the string representation of the current vectorial point.
     * @return the string representation of the current vectorial point.
     **/
    public String toString()
    {
        return "{ " + x + " , " + y + " }";
    }
}

Runge-Kutta à n dimensions :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
 
/*
 * RK4.java
 *
 * Created on 16 décembre 2006, 20:26
 *
 */
 
package org.umh.math.analysis.integrator;
 
import org.umh.math.linear.Vector;
 
/**
 * This class models a RK4 object that will be able to integrate vectorial
 * functions, i.e functions from R to Rn by the Runge-Kutta algorithm of 
 * order 4.
 * @author Absil Romain.
 */
public class RK4 implements RNIntegrator
{
 
    private RKFunction f;
    private double t0;
    private double t1;
    private Vector u0;
 
    /**
     * Construcs a new RK4 object with the specified function 
     * to integrate on the interval [t0, t1] and wich value on t0 is u0.
     * <p>i.e. Considering the following Cauchy problem : 
     * <p> d u(t) = f(t, u(t)) and f(t0) = u0     (with u : R -> Rn)
     * <p> the algotihm will return the function f integrated on the 
     * interval [t0, t1].
     * @param f the function to integrate.
     * @param t0 the lower born of the interval.
     * @param t1 the higher born of the interval.
     * @param u0 the function evaluated in t0.
     **/
    public RK4(RKFunction f, double t0, double t1, Vector u0)
    {
        this.f = f;
        this.t0 = t0;
        this.t1 = t1;
        this.u0 = u0;
    }
 
    /**
     * Returns the function integrated on n nodes in the interval [t0, t1].
     * The return type is some vectorial points of the integrated function.
     * @param n the number of nodes on wich you want the function to be 
     * integrated.
     * @return the function integrated on n nodes in the interval [t0, t1].
     **/
    public VectorialPoint[] integrate(int n)
    {
        double h = (t1 - t0) / (n - 1);
        VectorialPoint[] suite = new VectorialPoint[n];
        suite[0] = new VectorialPoint(t0, u0);
        Vector ui = u0;
 
        for(int i = 0 ; i < n ; i++)
        {
            double ti = t0 + (h * i);
            Vector k1 = f.eval(ti , ui).multByScal(h);
            Vector k2 = f.eval(ti + h/2. , 
                    ui.add(k1.multByScal(0.5))).multByScal(h);
            Vector k3 = f.eval(ti + h/2., 
                    ui.add(k2.multByScal(0.5))).multByScal(h);
            Vector k4 = f.eval(ti + h, ui.add(k3)).multByScal(h);
            ui = ui.add( 
                    (k1.add(k2.multByScal(2)).add(k3.multByScal(2)).add(k4))
                        .multByScal(1/6.) );
            suite[i] = new VectorialPoint(ti, ui);
        }
        return suite;
    }
 
}