API tools faq
paste
Ladies_Man

num meth sent

Ladies_Man
Sep 25th, 2016
104
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 9.01 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <math.h>
  3. #include "stdio.h"
  4.  
  5. using namespace std;
  6.  
  7. #define xtype float
  8. #define N 10
  9.  
  10. xtype _a = 2.477520;
  11. xtype _b = 0.232206;
  12.  
  13.  
  14. xtype xs[] = {0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5};
  15. xtype ys[] = {2.78, 3.13, 3.51, 3.94, 4.43, 4.97, 5.58, 6.27, 7.04, 7.91, 8.89};
  16.  
  17. 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};
  18. 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.39, 8.89};
  19.  
  20. xtype xsNew1[2*N+1];
  21. xtype ysNew1[2*N+1];
  22.  
  23. xtype f(xtype x) { return _a * exp(_b*x); }
  24.  
  25. xtype h = xs[1] - xs[0];
  26.  
  27. typedef struct spline {
  28.     xtype a[N];
  29.     xtype b[N];
  30.     xtype c[N];
  31.     xtype d[N];
  32. } _spline;
  33.  
  34. typedef struct bispline {
  35.     xtype x[N+3];
  36. } _bispline;
  37.  
  38. typedef struct akame {
  39.     xtype a[N+4];
  40.     xtype b[N+4];
  41.     xtype c[N+4];
  42.     xtype d[N+4];
  43. } _akame;
  44.  
  45. void solveMatrix(int n, xtype* a, xtype* b, xtype* c, xtype* d, xtype* x) {
  46.  
  47.     int i;
  48.     xtype alpha[n], beta[n];
  49.     alpha[0] = -c[0] / b[0];
  50.     beta[0] = d[0] / b[0];
  51.  
  52.     for (i = 1; i < n-1; ++i) {
  53.         alpha[i] = (-c[i]) / (a[i] * alpha[i-1] + b[i]);
  54.         beta[i] = (d[i] - a[i] * beta[i-1]) / (a[i] * alpha[i-1] + b[i]);
  55.     }
  56.  
  57.     beta[n-1] = (d[n-1]-a[n-1] * beta[n-2]) / (a[n-1] * alpha[n-2] + b[n-1]);
  58.     x[n-1] = beta[n-1];
  59.  
  60.     for (i = n-2; i >= 0; --i) {
  61.         x[i] = alpha[i] * x[i+1] + beta[i];
  62.     }
  63. }
  64.  
  65. void solveMatrixNew(int len, xtype* a, xtype* b, xtype* c, xtype* d, xtype* x) {
  66.  
  67.     int size = len;
  68.     int n = size - 1;
  69.     int i;
  70.  
  71.     xtype alpha[size], beta[size];
  72.  
  73.     alpha[0] = -c[0] / b[0];
  74.     beta[0] = d[0] / b[0];
  75.  
  76.     for (i = 1; i <= n; i++) {
  77.         alpha[i] = (-c[i]) / (a[i] * alpha[i-1] + b[i]);
  78.         beta[i] = (d[i] - a[i] * beta[i-1]) / (a[i] * alpha[i-1] + b[i]);
  79.     }
  80.  
  81.     //beta[n-1] = (d[n-1]-a[n-1] * beta[n-2]) / (a[n-1] * alpha[n-2] + b[n-1]);
  82.     x[n+1] = beta[n];
  83.  
  84.     for (i = n; i >= 1; i--) {
  85.         x[i] = alpha[i-1] * x[i+1] + beta[i-1];
  86.     }
  87.     x[0] = 2*x[1] - x[2];
  88.     x[n+2] = 2*x[n+1] - x[n];
  89. }
  90.  
  91.  
  92.  
  93.  
  94. spline populateSpline(int n) {
  95.  
  96.     xtype ak[n-3], bk[n-2], ck[n-3], dk[n-2], ek[n];
  97.     xtype a[n], b[n], c[n], d[n];
  98.  
  99.     xtype h = xs[1] - xs[0];
  100.  
  101.     for (int i = 0; i < n - 2; i++) {
  102.  
  103.         if (i < n-3)
  104.             ak[i] = ck[i] = 1;
  105.  
  106.         bk[i] = 4;
  107.         dk[i] = 3 / (h * h) * (ys[i+2] - 2*ys[i+1] + ys[i]);
  108.     }
  109.  
  110.     solveMatrix(n-2, ak, bk, ck, dk, ek);
  111.  
  112.  
  113.     c[0] = c[n-1] = 0;
  114.     for (int i = 1; i < n-1; i++)
  115.         c[i] = ek[i-1];
  116.  
  117.     for (int i = 0; i < n-1; i++) {
  118.         a[i] = ys[i];
  119.         b[i] = (ys[i+1] - ys[i]) / h - (h/3) * (c[i+1] + 2*c[i]);
  120.         d[i] = (c[i+1] - c[i]) / (3*h);
  121.     }
  122.  
  123.  
  124.     spline spZ;
  125.  
  126.     for (int i = 0; i < n; i++) {
  127.         spZ.a[i] = a[i];
  128.         spZ.b[i] = b[i];
  129.         spZ.c[i] = c[i];
  130.         spZ.d[i] = d[i];
  131.     }
  132.     return spZ;
  133. }
  134.  
  135. xtype countSplineAt(_spline spline, xtype x) {
  136.  
  137.     //S_j[x] = a_j + b_j(x-xj) + c_j(x-xj)^2 + d_j(x-xj)^3
  138.     int i;
  139.     for (i = 0; i < N; ++i)
  140.         if (x >= xs[i] && x <= xs[i+1])
  141.             break;
  142.  
  143.     if (N == i) i--;
  144.  
  145.     //printf("cSA: %.6f\n", xsNew[i]);
  146.     //printf("[%d] %.6f %.6f \n",i, x, xsNew[i]);
  147.  
  148.     return  spline.a[i] +
  149.             spline.b[i]*(x - xs[i]) +
  150.             spline.c[i]*(x - xs[i])*(x - xs[i]) +
  151.             spline.d[i]*(x - xs[i])*(x - xs[i])*(x - xs[i]);
  152. }
  153.  
  154.  
  155. xtype nFunc(int q, int k, xtype t) {
  156.  
  157.     if (1 == q) {
  158.         if (t >= xs[k] && t <= xs[k+1])
  159.             return 1;
  160.         else
  161.             return 0;
  162.     }
  163.  
  164.     xtype n1 = nFunc(q-1, k, t) * (t - xs[k]) / (xs[k+q-1] - xs[k]);
  165.     xtype n2 = nFunc(q-1, k+1, t) * (xs[k+q] - t) / (xs[k+q] - xs[k+1]);
  166.  
  167.     return n1 + n2;
  168. }
  169.  
  170.  
  171.  
  172. bispline populateBiSpline(int len) {
  173.  
  174.     int size = len;
  175.     int n = size - 1;
  176.  
  177.     xtype a[size], b[size], c[size], d[size], e[size+2];
  178.  
  179.     for (int i = 1; i < n; i++) {
  180.         a[i] = 1;
  181.         b[i] = 4;
  182.         c[i] = 1;
  183.         d[i] = 6*ys[i];
  184.     }
  185.  
  186.     a[0] = 0;
  187.     b[0] = 1;
  188.     c[0] = 0;
  189.     d[0] = ys[0];
  190.  
  191.     a[n] = 0;
  192.     b[n] = 1;
  193.     c[n] = 0;
  194.     d[n] = ys[n-1];
  195.  
  196.     for (int i = 0; i < n; i++) {
  197.         d[i] /= (h*h*h);
  198.     }
  199.  
  200.     solveMatrixNew(len, a, b, c, d, e);
  201.  
  202.  
  203.     _bispline bsp;
  204.  
  205.  
  206.     for (int i = 0; i < size+2; i++)
  207.         bsp.x[i] = e[i];
  208.  
  209.     return bsp;
  210. }
  211.  
  212. xtype baseX(int k, xtype x, xtype kh) {
  213.     if (0 == k) {
  214.  
  215.         xtype xz = xs[0];
  216.  
  217.         if (xz + 2*h <= x) return 0;
  218.         if (xz + h <= x && x < xz + 2*h) return pow(2*h + xz - kh - (x - kh), 3)/6;
  219.         if (xz <= x && x < xz + h)
  220.             return 2*pow(h, 3)/3 - pow(-xz + kh + (x - kh), 2)*(2*h + xz - kh - (x - kh))/2;
  221.         if (xz - h <= x && x < xz)
  222.             return 2*pow(h, 3)/3 - pow(-xz + kh + (x - kh), 2)*(2*h - xz + kh + (x - kh))/2;
  223.         if (xz - 2*h <= x && x < xz - h) return pow(2*h - xz + kh + (x - kh), 3)/6;
  224.         if (x < xz - 2*h) return 0;
  225.  
  226.     } else {
  227.         return baseX(0, x - k*h, -k*h);
  228.     }
  229.  
  230. }
  231.  
  232. xtype countBiSplineAt(_bispline bsp, xtype val) {
  233.  
  234.     int i;
  235.  
  236.     for (i = 0; i < N+2; ++i)
  237.         if (xs[i] <= val && val <= xs[i+1])
  238.             break;
  239.  
  240.     if (N+2 == i) i--;
  241.  
  242.     //printf("%.2f < %.2f < %.2f\n", xs[i], val, xs[i+1]);
  243.     //if (xs[i] <= val && val <= xs[i+1]) {
  244.         return  bsp.x[i + 0] * baseX(i - 1, val, 0) +
  245.                 bsp.x[i + 1] * baseX(i + 0, val, 0) +
  246.                 bsp.x[i + 2] * baseX(i + 1, val, 0) +
  247.                 bsp.x[i + 3] * baseX(i + 2, val, 0);
  248.     //} else {
  249.     //    return -1;
  250.     //}
  251. }
  252.  
  253.  
  254. akame populateAkame(int n) {
  255.  
  256.     int i;
  257.     xtype m[n + 4];
  258.     //-2;-1;    0.. n-1    ;n;n+1
  259.     for (i = 2; i < n + 2; i++) {
  260.         m[i] = (ys[i+1-2] - ys[i-2]) / (xs[i+1-2] - xs[i-2]);
  261.     }
  262.  
  263.     m[0] = 3*m[2] - 2*m[3]; //m[-2]=m[0]-m[1]
  264.     m[1] = 2*m[2] - m[3];
  265.     m[n+2] = 2*m[n+1] - m[n];
  266.     m[n+3] = 3*m[n];
  267.  
  268.     xtype ne, alpha_i, h_i;
  269.     xtype t_l[n], t_r[n];
  270.  
  271.     for (i = 2; i < n + 2; i++) {
  272.         ne = fabs(m[i+1] - m[i]) + fabs(m[i-1] - m[i-2]);
  273.  
  274.         if (ne > 0) {
  275.             alpha_i = fabs(m[i-1] - m[i-2]) / ne;
  276.             t_l[i] = m[i-1] + alpha_i*(m[i] - m[i-1]);
  277.             t_r[i] = t_l[i];
  278.         } else {
  279.             t_l[i] = m[i-1];
  280.             t_r[i] = m[i];
  281.         }
  282.     }
  283.  
  284.     akame ga_spline;
  285.     for (i = 2; i < n+2; i++) {
  286.         h_i = xs[i+1-2] - xs[i-2];
  287.  
  288.         ga_spline.a[i] = ys[i-2];
  289.         ga_spline.b[i] = t_r[i];
  290.         ga_spline.c[i] = (3*m[i] - 2*t_r[i] - t_l[i+1]) / h_i;
  291.         ga_spline.d[i] = (t_r[i] + t_l[i+1] - 2*m[i]) / (h_i*h_i);
  292.     }
  293.  
  294.     return ga_spline;
  295. }
  296.  
  297. xtype countAkameAt(akame ak, xtype x, int i) {
  298.  
  299.     for (i = 0; i <= N; ++i)
  300.         if (x >= xs[i] && x <= xs[i+1])
  301.             break;
  302.  
  303.     return  ak.a[i+2] +
  304.             ak.b[i+2]*(x - xs[i]) +
  305.             ak.c[i+2]*(x - xs[i])*(x - xs[i]) +
  306.             ak.d[i+2]*(x - xs[i])*(x - xs[i])*(x - xs[i]);
  307. }
  308.  
  309. int f1() {
  310.     char ch, ch1 = 'a', ch2 = 'b', ch3 = 'c';
  311.     ch = ch1 + ch2 + ch3;
  312.     //f1 = int(ch);
  313.     //return f1;
  314. }
  315.  
  316. int sum(int start, int end) {
  317.     if (end == start)
  318.         return 1;
  319.     else
  320.         return end + sum(start, end - 1);
  321. }
  322.  
  323. int main()
  324. {
  325.     for (int i = 0; i <= 2*N; i++) {
  326.         xtype val = (i % 2) == 0 ? xsNew[i] : xsNew[i-1] + 0.15;
  327.         xsNew1[i] = val;
  328.         printf("%.2f, ", val);
  329.     }
  330.     printf("\n===\n");
  331.     for (int i = 0; i <= 2*N; i++) {
  332.         xtype val = (i % 2) == 0 ? ysNew[i] : f(xsNew1[i]);
  333.         ysNew1[i]= val;
  334.         printf("%.2f, ", val);
  335.     }
  336.  
  337.  
  338.     printf("populating Spline\n");
  339.     spline sp = populateSpline(N+1);
  340.  
  341.     printf("populating Bispline\n");
  342.     bispline bs = populateBiSpline(N+2);
  343.  
  344.     printf("populating Akame ga Spline\n");
  345.     akame as = populateAkame(N);
  346.  
  347.     for ( int i = 1; i <= 2*N; i+=2) {
  348.  
  349.         xtype sVal = countSplineAt(sp, xsNew1[i]);
  350.         xtype bsVal = countBiSplineAt(bs, xsNew1[i]);
  351.         xtype asVal = countAkameAt(as, xsNew1[i], i);
  352.  
  353.         printf("x[%d]=%.2f y=%.5f \ts=%.5f d=%.5f \tbi=%.5f d=%.5f \ta=%.5f d=%.5f\n",
  354.                i, xsNew1[i], ysNew1[i],
  355.                sVal, fabs(ysNew1[i] - sVal),
  356.                bsVal, fabs(ysNew1[i] - bsVal),
  357.                asVal, fabs(ysNew1[i] - asVal)
  358.                );
  359.         //latex graph
  360.         /*printf("%.2f %.5f %.5f %.5f %.5f\n",
  361.                xsNew1[i], ysNew1[i],
  362.                sVal, //fabs(ysNew1[i] - sVal),
  363.                bsVal, //fabs(ysNew1[i] - bsVal),
  364.                asVal//, fabs(ysNew1[i] - asVal)
  365.                );*/
  366.         //latex
  367.         /*printf("%.2f & %.6f & %.6f & %.6f & %.6f & %.6f & %.6f & %.6f \\\\\n",
  368.                 xsNew1[i], ysNew1[i],
  369.                 sVal, fabs(ysNew1[i] - sVal),
  370.                 bsVal, fabs(ysNew1[i] - bsVal),
  371.                 asVal, fabs(ysNew1[i] - asVal)
  372.                 );
  373.         printf("\\hline\n");*/
  374.     }
  375.     return 0;
  376. }
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!