Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

Untitled

By: a guest on May 2nd, 2012  |  syntax: None  |  size: 4.26 KB  |  hits: 17  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. /*  MATH1155 C programming assignment */
  2. /*  file: intersect.c                 */
  3. /*  created by: Ian Grundy 01/04/12   */
  4. /*  modified by:                      */
  5.  
  6. # include <stdio.h>
  7. # include <stdlib.h>
  8. # include <math.h>
  9.  
  10. # define PI 3.1415926535897932385
  11.  
  12. /*  Function prototypes */
  13.  
  14. double f(double x);
  15. double dfdx(double x);
  16.  
  17. double si(double x);
  18. double diffsi(double x);
  19.  
  20. double myrandom();
  21. double heuristic(double left, double right, int maxiter);
  22. double bisection(double left, double right, double epsilon);
  23. double secant(double left, double right, double epsilon, int maxiter);
  24. double newton(double initial, double epsilon, int maxiter);
  25.  
  26. int main(void)
  27. {
  28.  
  29.    double eps = 1.0e-12;
  30.    int maxiter;
  31.    double b, n, s, h;
  32.  
  33.    maxiter = 5000;
  34.    h = heuristic(0, 3.0, maxiter);
  35.  
  36.    maxiter = 500;
  37.    b = bisection(0, 3.0, eps);
  38.    s = secant(0, 3.0, eps, maxiter);
  39.    n = newton(1.5, eps, maxiter);
  40.  
  41.    printf ("   heuristic       bisection       secant       newton\n");
  42.    printf ("%15.10f %15.10f %15.10f %15.10f %f\n", h, b, s, n, si(1.5));
  43.  
  44.    return 0;
  45.  
  46. }
  47.  
  48. /*  Function Definitions */
  49.  
  50.     double factorial(int n)
  51.     {
  52.         int i;
  53.         double product = 1;
  54.         for ( i=2; i <= n; i=i+1) product = product * i;
  55.         return product;
  56.     }
  57.  
  58.  
  59.     double si(double x)
  60.     {
  61.         int i;
  62.         double z;
  63.         double sum=0;
  64.         for(i=1;;i++)
  65.         {
  66.             z = ((pow((-1),i))*(pow(x,(2*i+1))))/((2*i+1)*(factorial(2*i+1)));
  67.             if (fabs(z) < exp(-12))
  68.             {
  69.                 break;
  70.             }
  71.             else
  72.             {
  73.                 sum = sum + z;
  74.             }
  75.  
  76.         }
  77.         return sum;
  78.     }
  79.  
  80.     double diffsi(double x)
  81.     {
  82.         double sum = 1;
  83.         int i;
  84.         double z;
  85.  
  86.         for(i=0;;i++)
  87.         {
  88.             z = ((pow((-1),i))*(pow(x,(2*i+2))))/(factorial(2*i+1));
  89.             if (fabs(z) < exp(-12))
  90.             {
  91.                 break;
  92.             }
  93.             else
  94.             {
  95.                 sum = sum + z;
  96.             }
  97.         }
  98.         return sum;
  99.     }
  100.  
  101.  
  102.  
  103.    double f(double x)
  104.    {
  105.       return si(x)-PI/2.0;
  106.    }
  107.  
  108.    double dfdx(double x)
  109.    {
  110.       return diffsi(x);
  111.    }
  112.  
  113.  
  114.    double myrandom()
  115.    {
  116.       return (double) rand() / (double) RAND_MAX;
  117.    }
  118.  
  119.    double heuristic(double left, double right, int maxiter)
  120.    {
  121.  
  122.       /* simulated annealing function minimiser */
  123.  
  124.       double rad = (right - left)/100.0;
  125.       double r = 0.99;
  126.       double T = 100;
  127.       double prob;
  128.  
  129.       int iter = 0;
  130.       double bestx, currentx, newx;
  131.       double besty, currenty, newy;
  132.  
  133.       currentx = (left + right)/2.0;
  134.       currenty = fabs(f(currentx));
  135.  
  136.       bestx = currentx;
  137.       besty = currenty;
  138.  
  139.       for (iter = 1; iter < maxiter+1; iter++)
  140.       {
  141.           newx = currentx + (2*myrandom()-1)*rad;
  142.           newy = fabs(f(newx));
  143.  
  144.           if ( newy < currenty )
  145.           {
  146.              /* accept if better */
  147.              currentx = newx;
  148.              currenty = newy;
  149.           }
  150.           else
  151.           {
  152.              /* possibly accept if worse */
  153.              prob = exp(-(newy-currenty)/T);
  154.              if ( myrandom() < prob)
  155.              {
  156.                 currentx = newx;
  157.                 currenty = newy;
  158.              }
  159.           }
  160.  
  161.           if ( currenty < besty )
  162.           {
  163.              /* store absolute best found */
  164.              bestx = currentx;
  165.              besty = currenty;
  166.           }
  167.  
  168.           T = r*T;
  169.  
  170.       }
  171.  
  172.       return bestx;
  173.    }
  174.  
  175.    double bisection(double left, double right, double epsilon)
  176.    {
  177.      double midpoint, fmid, fleft, fright;
  178.  
  179.      fleft = f(left);
  180.      fright = f(right);
  181.  
  182.      if (fleft*fright > 0) return -99;
  183.  
  184.      while (fabs(right - left) > 2*epsilon)
  185.      {
  186.            midpoint = (left + right)/2.0;
  187.            fmid = f(midpoint);
  188.  
  189.            if (fmid*fleft >= 0) left = midpoint;
  190.            if (fmid*fright > 0) right = midpoint;
  191.      }
  192.  
  193.      return (left + right)/2.0;
  194.    }
  195.  
  196.    double secant(double left, double right, double epsilon, int maxiter)
  197.    {
  198.      return (left + right)/2.0;
  199.    }
  200.  
  201.    double newton(double initial, double epsilon, int maxiter)
  202.    {
  203.       return initial;
  204.    }
create a new version of this paste RAW Paste Data