Opdrachten numerieke methoden, serie 3

(1)

Opdrachten numerieke methoden, serie 3

Opdracht 1: De Nulpunten van Legendre Polynomen

p(x) = c

N

Y

n=1

(x + a

n

)

p

⁰

(x) = c

N

X

n=1 N

Y

m=1,m6=n

(x + a

m

)

p

⁰⁰

(x) = c

N

X

n=1 N

X

m=1,m6=n N

Y

p=1,p6=m,p6=n

(x + a

p

)

Als x > x

^∗

, dan is Q

N

n=1

(x + a

_n

) > 0, daar a

_n

de nulpunten van het polynoom zijn en x dus groter is dan het grootste nulpunt. Dan

p(x) · p

⁰⁰

(x) = c

²

N

X

n=1 N

X

m=1,m6=n N

Y

p=1,p6=n,p6=m

(x + a

p

)

N

Y

q=1

(x + a

q

).

Zowel c

²

als (x + a

i

) zijn groter dan nul, dus p(x) · p

⁰⁰

(x) > 0.

Het programma

/*

* legendre.c

*

* by Willem van Engen and Hjalmar Mulders

*/

#include <stdio.h>

#include <math.h>

#include <err.h>

typedef double (*t iterate function)(const double,double*,const double*,const int); 10

#define N 7

#define POLYARRAY {0.,−35./16.,0.,315./16.,0.,−693./16.,0.,429./16.}

/*

* double horner(const double x, double *b, const double *a, const int n)

*

* Description:

* Computes the horner coefficients of the given polynomial in point x and computes p(x)

* 20

* Arguments:

* const double x - point of evaluation

* double b - destination array of Horner coefficients [0. .n]

* const double a - source array of coefficients of p(x) [0. .n]

* const int n - degree of polynomial

* Returns:

* the evaluated polynomial in point x

*/

double

horner(const double x, double *b, const double *a, const int n){ ³⁰

int i;

b[n]=a[n];

(2)

for (i=n−1;i>=0;i−−) { b[i]=b[i+1]*x+a[i];

}

return b[0];

}

/* 40

* double newtonstep(const double x, double *b, const double *a, const int n)

*

* Description:

* Compute the derivative of the polynomial in point x.

*

* Arguments:

* const double x - point of evaluation

* double b - dest. array of derivative’s coefficients [0. .n]

* const double a - source array of polynomial coefficients [0. .n]

* const int n - degree of polynomial 50

* Returns:

* The value of the derivative in point x

*/

double

newtonstep(const double x, double *b, const double *a, const int n){ double peval;

double c[N+1];

peval=horner(x,b,a,n);

horner(x,c,b,n); 60

return x− peval/c[1];

}

/*

* iterate fixed point(double *x, int *K, double *b, const double *a, const int n, const double eps, const int m, t iterate function)

*

* Description:

* Iterate to the largest root of the given polynomial

* 70

* Arguments:

* double *x - array with startvalue x[0] and destination for the computed approximations

* int *K - number of values in *x

* double *b - destination horner coefficients of the polynomial [0. .n1]

* const double *a - array with coefficients of the polynomial

* const int n - degree of polynomial

* const double eps - stop when this precision is reached

* const int m - stop when done this many iterations

* t iterate fuction g - fixed point function

* 80

* Returns:

* Last computed value for the largest root

*/

double

iterate fixed point(double *x, int *k, double *b, const double *a, const int n, const double eps, const int m, t iterate function g) {

for ((*k)=1; (m−1)>=(*k); (*k)++) { x[*k]=g(x[(*k)−1],b,a,n);

if (fabs(x[*k]−x[(*k)−1])<eps) break;

} ⁹⁰

return x[*k>=m?m−1:*k];

}

int main(int argc, char *argv[ ]){ int m=100;

double eps=1e−12;

double A[N+1]=POLYARRAY;

double B[N+1];

double *a=A,*b=B,*c; 100

(3)

int K,i,j;

double *x;

if ( (x=(double*)malloc(sizeof (double)*m))==NULL) return 1;

x[0]=1.5;

for (i=0;i<N/2;i++){

printf("%16.16e:\n---\n", iterate fixed point(x,&K,b,a,N−i,eps,m,newtonstep) );

for (j=0;j<K;j++){

printf(" x[%2.2d] = %16.16e\n",j,x[j]); 110

}

printf("\n\n");

c=a; a=&(b[1]); b=c;

}

free(x);

return 0;

}

De output van dit programma is:

9.4910791234275860e-01:

---

x[00] = 1.5000000000000000e+00 x[01] = 1.3358889406323233e+00 x[02] = 1.2040473175706508e+00 x[03] = 1.1016671338433446e+00 x[04] = 1.0268724008367545e+00 x[05] = 9.7860406647878473e-01 x[06] = 9.5519091852670202e-01 x[07] = 9.4943682214155622e-01 x[08] = 9.4910894317751215e-01 x[09] = 9.4910791235292558e-01 x[10] = 9.4910791234275849e-01 7.4153118559939468e-01:

---

x[00] = 1.5000000000000000e+00 x[01] = 1.2662566736985428e+00 x[02] = 1.0812505162996384e+00 x[03] = 9.3966043078549399e-01 x[04] = 8.3816034767385239e-01 x[05] = 7.7507468402843216e-01 x[06] = 7.4721699085777893e-01 x[07] = 7.4173070605518476e-01 x[08] = 7.4153144252668113e-01 x[09] = 7.4153118559982145e-01 4.0584515137561883e-01:

---

x[00] = 1.5000000000000000e+00 x[01] = 1.1621337678560417e+00 x[02] = 9.0031622631797603e-01 x[03] = 7.0223621044698903e-01 x[04] = 5.5935056511268888e-01 x[05] = 4.6649930368820031e-01 x[06] = 4.1957536455921091e-01 x[07] = 4.0675794491967571e-01 x[08] = 4.0584954394903450e-01 x[09] = 4.0584515147798517e-01 x[10] = 4.0584515137561883e-01

Deze waarden komen overeen met de wortels uit het dictaat op pagina 36.

Opdracht 2: Lineaire stelsels

Voor het beschouwde circuit wordt op de knooppunten de som der stromen nul gesteld. Hier wordt een voorbeeld van gegeven op v

1

:

10 − v

1

R

1

− v

₁

− v

4

R

4

− v

₁

− v

3

R

3

= 0.

(4)

Dit vereenvoudigd levert:

3v

1

+ v

3

+ v

4

= 10

Uiteindelijk leveren de vergelijkingen voor deze knooppunten het volgende stelsel op:

3v

1

−v

3

−v

4

= 10

3

2

v

2

−

¹₂

v

4

−

¹₂

v

5

= 0

−v

1 3

2

v

3

−

¹₂

v

6

= 0

−v

1

−

¹₂

v

2

3v

4

−

¹₂

v

6

−v

7

= 0

−

¹₂

v

2 3

2

v

5

−v

7

= 0

−

¹₂

v

3

−

¹₂

v

4 3

2

v

6

= 0

−v

4

−v

5

3v

7

= 0

Dat een matrix diagonaal is dominant, wordt bewezen met

n

X

i=1,i6=j

|a

ij

| ≥ |a

jj

| j = 1, . . . , n

Het is gemakkelijk na te gaan dat dit voor deze matrix zo is; de diagonaalcel in iedere kolom is groter dan de som van de rest in de kolom.

Het programma

/*

* gauss.c

*

* by Willem van Engen and Hjalmar Mulders

*/

#include <stdio.h>

#define N 7

#define MATRIXA { \

{ 3, 0, −1, −1, 0, 0, 0}, \ ¹⁰

{ 0, 1.5, 0,−0.5,−0.5, 0, 0}, \ { −1, 0, 1.5, 0, 0,−0.5, 0}, \ { −1,−0.5, 0, 3, 0,−0.5, −1}, \ { 0,−0.5, 0, 0, 1.5, 0, −1}, \ { 0, 0,−0.5,−0.5, 0, 1.5, 0}, \ { 0, 0, 0, −1, −1, 0, 3} }

#define VECTORB{10, 5, 0, 0, 0, 0, 0 }

#define TRUE 0 20

#define FALSE 1 typedef char bool;

/*

* void gauss(int n, double **A, double *b, double *x, bool acc)

*

* Description:

* Gauss-elimination of Ax=b.

*

* Arguments: 30

* int n - dimension of matrix A

* double **a - the source matrix A

* double *b - the source vector b (size=n)

* Returns:

* TRUE if succesful, FALSE if failed

*/

bool

gauss(int n, double A[N][N], double b[N], double x[N]){ int i,j,k;

double l; 40

(5)

for (k=0;k<n;k++){ for (i=k+1;i<n;i++){

l=A[i][k]/A[k][k];

for (j=k+1;j<n;j++){ A[i][j]−=l*A[k][j];

}

b[i]−=l*b[k];

}

} ⁵⁰

for (k=n−1;k>=0;k−−) { for (j=k+1;j<n;j++){ b[k]−=x[j]*A[k][j];

}

x[k]=b[k]/A[k][k];

} }

int 60

main(int argc, char *argv[ ]){ double a[N][N] = MATRIXA;

double b[N] = VECTORB;

double x[N];

int i;

gauss(N,a,b,x);

for (i=0;i<N;i++) printf(" x[%1.1d] = %16.16f\n",i,x[i]);

return 0;

} ⁷⁰

Dit programma wordt getest met het stelsel voor het electrische circuit. De output van het programma is dan:

x[0] = 6.9444444444444438 x[1] = 6.3888888888888893 x[2] = 5.8333333333333330 x[3] = 5.0000000000000009 x[4] = 4.1666666666666679 x[5] = 3.6111111111111107 x[6] = 3.0555555555555562

Dit komt overeen met een berekening in matlab.

Willem van Engen (wvengen@stack.nl) and Hjalmar Mulders (H.C.J.Mulders@student.tue.nl)

Homepage: http://willem.n3.net