Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <math.h>
- #include "stdio.h"
- using namespace std;
- #define xtype double
- #define N 10
- xtype _a = 2.477520;
- xtype _b = 0.232206;
- xtype xs[] = {0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5};
- xtype ys[] = {2.78, 3.13, 3.51, 3.94, 4.43, 4.97, 5.58, 6.27, 7.04, 7.91, 8.89};
- xtype xsNew[] = {0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00, 2.25, 2.50, 2.75, 3.00, 3.25, 3.50, 3.75, 4.00, 4.25, 4.50, 4.75, 5.00, 5.25, 5.50};
- xtype ysNew[] = {2.78, 2.95, 3.13, 3.31, 3.51, 3.72, 3.94, 4.18, 4.43, 4.69, 4.97, 5.27, 5.58, 5.92, 6.27, 6.65, 7.04, 7.47, 7.91, 8.38, 8.89};
- xtype f(xtype x) { return _a * exp(_b*x); }
- xtype h = xs[1] - xs[0];
- typedef struct spline {
- xtype a[N];
- xtype b[N];
- xtype c[N];
- xtype d[N];
- } _spline;
- typedef struct bispline {
- xtype x[N+3];
- } _bispline;
- void solveMatrix(int n, xtype* a, xtype* b, xtype* c, xtype* d, xtype* x) {
- int i;
- xtype alpha[n], beta[n];
- alpha[0] = -c[0] / b[0];
- beta[0] = d[0] / b[0];
- for (i = 1; i < n-1; ++i) {
- alpha[i] = (-c[i]) / (a[i] * alpha[i-1] + b[i]);
- beta[i] = (d[i] - a[i] * beta[i-1]) / (a[i] * alpha[i-1] + b[i]);
- }
- beta[n-1] = (d[n-1]-a[n-1] * beta[n-2]) / (a[n-1] * alpha[n-2] + b[n-1]);
- x[n-1] = beta[n-1];
- for (i = n-2; i >= 0; --i) {
- x[i] = alpha[i] * x[i+1] + beta[i];
- }
- }
- void solveMatrixNew(int len, xtype* a, xtype* b, xtype* c, xtype* d, xtype* x) {
- int size = len;
- int n = size - 1;
- int i;
- xtype alpha[size], beta[size];
- alpha[0] = -c[0] / b[0];
- beta[0] = d[0] / b[0];
- for (i = 1; i <= n; i++) {
- alpha[i] = (-c[i]) / (a[i] * alpha[i-1] + b[i]);
- beta[i] = (d[i] - a[i] * beta[i-1]) / (a[i] * alpha[i-1] + b[i]);
- }
- //beta[n-1] = (d[n-1]-a[n-1] * beta[n-2]) / (a[n-1] * alpha[n-2] + b[n-1]);
- x[n+1] = beta[n];
- for (i = n; i >= 1; i--) {
- x[i] = alpha[i-1] * x[i+1] + beta[i-1];
- }
- x[0] = 2*x[1] - x[2];
- x[n+2] = 2*x[n+1] - x[n];
- }
- _spline populateSpline(int n) {
- xtype ak[n-3], bk[n-2], ck[n-3], dk[n-2], ek[n];
- xtype a[n], b[n], c[n], d[n];
- xtype h = xs[1] - xs[0];
- for (int i = 0; i < n - 2; i++) {
- if (i < n-3)
- ak[i] = ck[i] = 1;
- bk[i] = 4;
- dk[i] = 3 / (h * h) * (ys[i+2] - 2*ys[i+1] + ys[i]);
- }
- solveMatrix(n-2, ak, bk, ck, dk, ek);
- c[0] = c[n-1] = 0;
- for (int i = 1; i < n-1; i++)
- c[i] = ek[i-1];
- for (int i = 0; i < n-1; i++) {
- a[i] = ys[i];
- b[i] = (ys[i+1] - ys[i]) / h - (h/3) * (c[i+1] + 2*c[i]);
- d[i] = (c[i+1] - c[i]) / (3*h);
- }
- _spline spZ;
- for (int i = 0; i < n; i++) {
- spZ.a[i] = a[i];
- spZ.b[i] = b[i];
- spZ.c[i] = c[i];
- spZ.d[i] = d[i];
- }
- return spZ;
- }
- xtype countSplineAt(_spline spline, xtype x) {
- //S_j[x] = a_j + b_j(x-xj) + c_j(x-xj)^2 + d_j(x-xj)^3
- int i;
- for (i = 0; i < N; ++i)
- if (x >= xs[i] && x <= xs[i+1])
- break;
- if (N == i) i--;
- //printf("cSA: %.6f\n", xsNew[i]);
- //printf("[%d] %.6f %.6f \n",i, x, xsNew[i]);
- return spline.a[i] +
- spline.b[i]*(x - xs[i]) +
- spline.c[i]*(x - xs[i])*(x - xs[i]) +
- spline.d[i]*(x - xs[i])*(x - xs[i])*(x - xs[i]);
- }
- xtype nFunc(int q, int k, xtype t) {
- if (1 == q) {
- if (t >= xs[k] && t <= xs[k+1])
- return 1;
- else
- return 0;
- }
- xtype n1 = nFunc(q-1, k, t) * (t - xs[k]) / (xs[k+q-1] - xs[k]);
- xtype n2 = nFunc(q-1, k+1, t) * (xs[k+q] - t) / (xs[k+q] - xs[k+1]);
- return n1 + n2;
- }
- bispline populateBiSpline(int len) {
- int size = len;
- int n = size - 1;
- xtype a[size], b[size], c[size], d[size], e[size+2];
- for (int i = 1; i < n; i++) {
- a[i] = 1;
- b[i] = 4;
- c[i] = 1;
- d[i] = 6*ys[i];
- }
- a[0] = 0;
- b[0] = 1;
- c[0] = 0;
- d[0] = ys[0];
- a[n] = 0;
- b[n] = 1;
- c[n] = 0;
- d[n] = ys[n];
- for (int i = 0; i < n; i++) {
- d[i] /= (h*h*h);
- }
- solveMatrixNew(len, a, b, c, d, e);
- _bispline bsp;
- for (int i = 0; i < size+2; i++)
- bsp.x[i] = e[i];
- return bsp;
- }
- xtype baseX(int k, xtype x, xtype kh) {
- if (0 == k) {
- xtype xz = xs[0];
- if (xz + 2*h <= x) return 0;
- if (xz + h <= x && x < xz + 2*h) return pow(2*h + xz - kh - (x - kh), 3)/6;
- if (xz <= x && x < xz + h)
- return 2*pow(h, 3)/3 - pow(-xz + kh + (x - kh), 2)*(2*h + xz - kh - (x - kh))/2;
- if (xz - h <= x && x < xz)
- return 2*pow(h, 3)/3 - pow(-xz + kh + (x - kh), 2)*(2*h - xz + kh + (x - kh))/2;
- if (xz - 2*h <= x && x < xz - h) return pow(2*h - xz + kh + (x - kh), 3)/6;
- if (x < xz - 2*h) return 0;
- } else {
- return baseX(0, x - k*h, -k*h);
- }
- }
- xtype countBiSplineAt(_bispline bsp, xtype val) {
- int i;
- for (i = 0; i < N+2; ++i)
- if (xs[i] <= val && val <= xs[i+1])
- break;
- if (N+2 == i) i--;
- //printf("%.2f < %.2f < %.2f\n", xs[i], val, xs[i+1]);
- //if (xs[i] <= val && val <= xs[i+1]) {
- return bsp.x[i + 0] * baseX(i - 1, val, 0) +
- bsp.x[i + 1] * baseX(i + 0, val, 0) +
- bsp.x[i + 2] * baseX(i + 1, val, 0) +
- bsp.x[i + 3] * baseX(i + 2, val, 0);
- /*} else {
- return -1;
- }*/
- }
- int main()
- {
- printf("populating spline\n");
- _spline sp = populateSpline(N+1);
- printf("populating bispline\n");
- _bispline bs = populateBiSpline(N+2);
- /*for (int i = 0; i <= 2*N; i++) {
- xtype val = xsNew[i];
- xtype bsVal = countBiSplineAt(bs, val);
- printf("[%d] x[i]:%.5f y:%.5f bs:%.5f :%.5f\n",
- i, val, ysNew[i], bsVal, fabs(ysNew[i] - bsVal));
- }*/
- for ( int i = 0; i <= 2*N; i++) {
- xtype sVal = countSplineAt(sp, xsNew[i]);
- xtype bsVal = countBiSplineAt(bs, xsNew[i]);
- printf("x[%d]=%.2f y=%.5f s=%.5f d=%.5f b=%.5f d=%.5f\n",
- i, xsNew[i], ysNew[i], sVal, fabs(ysNew[i] - sVal), bsVal, fabs(ysNew[i] - bsVal));
- //printf("x[%d]=%.2f y[%d]=%.6f bs[%d]=%.6f\n",
- // i, xsNew[i], i, ysNew[i], i, bsVal);
- }
- /* xtype s = 0;
- int k = 3;
- int q = 4;
- for ( int k = 0; k < N; k++) {
- s += nFunc(q, k, 3.0) * ys[k];
- }
- printf("s=%.6f", s);*/
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement