pilotpirx/numerik1

Aufgabe 2 fertig

aseyboldt committed Apr 26, 2012
1 parent d173540 commit 87eba82ccd05c9f08dbfaa28ba56cd0c54c9483c
Showing with 21 additions and 8 deletions.
1. +8 −2 aufgabe02.py
2. +11 −4 poly.c
3. +2 −2 poly.h
 @@ -8,18 +8,19 @@ _libpoly = np.ctypeslib.load_library('libpoly', '.') _libpoly.poly_de_casteljau_array.argtypes = [_ndpointer_double, + ctypes.c_int, ctypes.c_int, _ndpointer_double, _ndpointer_double, ctypes.c_int] -def poly_de_casteljau(points, t): +def poly_de_casteljau(points, t, n_deriv=0): points = np.asarray(points, dtype=np.double) t = np.asarray(t, dtype=np.double) out_all_dims = [] for b in points: out = np.empty(t.size, np.double) - _libpoly.poly_de_casteljau_array(b, len(b), t.ravel(), out, t.size) + _libpoly.poly_de_casteljau_array(b, len(b), n_deriv, t.ravel(), out, t.size) out_all_dims.append(out.reshape(t.shape)) return np.array(out_all_dims) @@ -33,6 +34,11 @@ def main(): plt.plot(*poly_de_casteljau([x[1:], y[1:]], t), alpha=0.5, linewidth=0.5) plt.plot(x, y, '--', alpha=0.5, linewidth=0.5) plt.plot(x, y, 'o', color='red') + + #dp_x, dp_y = poly_de_casteljau([x, y], t, 1) + #d2p_x, d2p_y = poly_de_casteljau([x, y], t, 2) + #plt.plot(t, dp_y) + #plt.plot(t, d2p_y) plt.show() if __name__ == '__main__':
15 poly.c
 @@ -32,7 +32,7 @@ void poly_neville_array(double *x, double *f, int n, double *t, double *out, int } -double poly_de_casteljau(double *b, double t, int n) +double poly_de_casteljau(double *b, double t, int n, int n_deriv) { double *array = malloc(n * sizeof(double)); if (!array) { @@ -46,23 +46,30 @@ double poly_de_casteljau(double *b, double t, int n) array[i] = b[i]; } - for (i = 1; i < n; i++) { + for (i = 1; i < n - n_deriv; i++) { for (j = 0; j < n - i; j++) { array[j] = t * array[j + 1] + (1 - t) * array[j]; } } + for (i = n - n_deriv; i < n; i++) { + for (j = 0; j < n - i; j++) { + array[j] = array[j + 1] - array[j]; + array[j] *= i; + } + } + tmp = array[0]; free(array); return tmp; } -void poly_de_casteljau_array(double *b, int n, double *t, double *out, int m) +void poly_de_casteljau_array(double *b, int n, int n_deriv, double *t, double *out, int m) { int i; for (i = 0; i < m; i++) { - out[i] = poly_de_casteljau(b, t[i], n); + out[i] = poly_de_casteljau(b, t[i], n, n_deriv); } }
 @@ -4,8 +4,8 @@ double poly_neville(double* x, double* f, double t, int n); void poly_neville_array(double *x, double *f, int n, double *t, double *out, int m); -double poly_de_casteljau(double *b, double t, int n); -void poly_de_casteljau_array(double *b, int n, double *t, double *out, int m); +double poly_de_casteljau(double *b, double t, int n, int n_deriv); +void poly_de_casteljau_array(double *b, int n, int n_deriv, double *t, double *out, int m); #endif