# Untitled

By: a guest on May 2nd, 2012  |  syntax: None  |  size: 4.26 KB  |  hits: 17  |  expires: Never
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.    }