Subversion Repositories programming

Rev

Rev 170 | Rev 173 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 170 Rev 171
Line 93... Line 93...
93 
{
 93 
{
94 
    float R[n+1][n+1];
 94 
    float R[n+1][n+1];
95 
    int i, j, k;
 95 
    int i, j, k;
96 
    float h, sum;
 96 
    float h, sum;
97 
 
 97 
 
- 
 
 98 
    /* Check for values that make us run wild, and refuse to run */
- 
 
 99 
    if (n > 25)
- 
 
 100 
    {
- 
 
 101 
        printf ("Cowardly refusing to run... You shouldn't need to use a\n"
- 
 
 102 
                "value of n > 25. Remember, the inner loop runs (2^n) - 1.\n");
- 
 
 103 
        return 0.0;
- 
 
 104 
    }
- 
 
 105 
 
98 
    h = b - a;
 106 
    h = b - a;
99 
    R[0][0] = (h/2.0) * (f(a) + f(b));
 107 
    R[0][0] = (h/2.0) * (f(a) + f(b));
100 
 
 108 
 
101 
    for (i=1; i<=n; i++)
 109 
    for (i=1; i<=n; i++)
102 
    {
 110 
    {
Line 156... Line 164...
156 
    float h, comp, abs, rel;
 164 
    float h, comp, abs, rel;
157 
    int i, n;
 165 
    int i, n;
158 
 
 166 
 
159 
    /* Euler Method Header */
 167 
    /* Euler Method Header */
160 
    printf ("Euler's Method\n");
 168 
    printf ("Euler's Method\n");
161 
    printf ("n       h               Computed        Abs Error       Rel Error\n");
 169 
    printf ("i       h               Computed        Abs Error       Rel Error\n");
162 
    printf ("====================================================================\n");
 170 
    printf ("====================================================================\n");
163
171 
 
164 
    /* Euler Method Table */
 172 
    /* Euler Method Table */
165 
    for (n=1; n<=25; n++)
 173 
    for (i=1; i<=25; i++)
166 
    {
 174 
    {
- 
 
 175 
        h = pow(2.0, -i);
167 
        h = (b - a) / n;
 176 
        n = (int)((b - a) / h);
168 
        comp = eval_euler (a, b, n, &func);
 177 
        comp = eval_euler (a, b, n, &func);
169 
        abs = fabs(real - comp);
 178 
        abs = fabs(real - comp);
170 
        rel = fabs(real - comp) / fabs(real);
 179 
        rel = fabs(real - comp) / fabs(real);
171 
 
 180 
 
172 
        printf ("%.2d\t%e\t%e\t%e\t%e\n", n, h, comp, abs, rel);
 181 
        printf ("%.2d\t%e\t%e\t%e\t%e\n", i, h, comp, abs, rel);
173 
    }
 182 
    }
174 
 
 183 
 
175 
    /* Trapezoid Method Header */
 184 
    /* Trapezoid Method Header */
176 
    printf ("\n\nTrapezoid Method\n");
 185 
    printf ("\n\nTrapezoid Method\n");
177 
    printf ("n       h               Computed        Abs Error       Rel Error\n");
 186 
    printf ("i       h               Computed        Abs Error       Rel Error\n");
178 
    printf ("====================================================================\n");
 187 
    printf ("====================================================================\n");
179
188 
 
180 
    /* Trapezoid Method Table */
 189 
    /* Trapezoid Method Table */
181 
    for (n=1; n<=25; n++)
 190 
    for (i=1; i<=25; i++)
182 
    {
 191 
    {
- 
 
 192 
        h = pow(2.0, -i);
183 
        h = (b - a) / n;
 193 
        n = (int)((b - a) / h);
184 
        comp = eval_trapezoid (a, b, n, &func);
 194 
        comp = eval_trapezoid (a, b, n, &func);
185 
        abs = fabs(real - comp);
 195 
        abs = fabs(real - comp);
186 
        rel = fabs(real - comp) / fabs(real);
 196 
        rel = fabs(real - comp) / fabs(real);
187 
 
 197 
 
188 
        printf ("%.2d\t%e\t%e\t%e\t%e\n", n, h, comp, abs, rel);
 198 
        printf ("%.2d\t%e\t%e\t%e\t%e\n", i, h, comp, abs, rel);
189 
    }
 199 
    }
190 
 
 200 
 
191 
    /* Romberg Method Header */
 201 
    /* Romberg Method Header */
192 
    printf ("\n\nRomberg's Method\n");
 202 
    printf ("\n\nRomberg's Method\n");
193 
    printf ("n       h               Computed        Abs Error       Rel Error\n");
 203 
    printf ("i       h               Computed        Abs Error       Rel Error\n");
194 
    printf ("====================================================================\n");
 204 
    printf ("====================================================================\n");
195
205 
 
196 
    /* Romberg Method Table */
 206 
    /* Romberg Method Table
- 
 
 207 
     *
- 
 
 208 
     * NOTE: Don't go higher than 3 for this one, since the Romberg method's inner
- 
 
 209 
     *       loop runs (2^24) - 1 = 16777216 times at i=3 and
- 
 
 210 
     *       (2^48) - 1 = 2.81474 * 10^14 times at i = 4.
- 
 
 211 
     */
197 
    for (n=1; n<=25; n++)
 212 
    for (i=1; i<=3; i++)
198 
    {
 213 
    {
- 
 
 214 
        h = pow(2.0, -i);
199 
        h = (b - a) / n;
 215 
        n = (int)((b - a) / h);
200 
        comp = eval_romberg (a, b, n, &func);
 216 
        comp = eval_romberg (a, b, n, &func);
201 
        abs = fabs(real - comp);
 217 
        abs = fabs(real - comp);
202 
        rel = fabs(real - comp) / fabs(real);
 218 
        rel = fabs(real - comp) / fabs(real);
203 
 
 219 
 
204 
        printf ("%.2d\t%e\t%e\t%e\t%e\n", n, h, comp, abs, rel);
 220 
        printf ("%.2d\t%e\t%e\t%e\t%e\n", i, h, comp, abs, rel);
205 
    }
 221 
    }
206 
 
 222 
 
207 
    /* Simpson Method Header */
 223 
    /* Simpson Method Header */
208 
    printf ("\n\nSimpson's Method\n");
 224 
    printf ("\n\nSimpson's Method\n");
209 
    printf ("n       h               Computed        Abs Error       Rel Error\n");
 225 
    printf ("i       h               Computed        Abs Error       Rel Error\n");
210 
    printf ("====================================================================\n");
 226 
    printf ("====================================================================\n");
211
227 
 
212 
    /* Simpson Method Table */
 228 
    /* Simpson Method Table */
213 
    for (n=1; n<=25; n++)
 229 
    for (i=1; i<=25; i++)
214 
    {
 230 
    {
- 
 
 231 
        h = pow(2.0, -i);
215 
        h = (b - a) / n;
 232 
        n = (int)((b - a) / h);
216 
        comp = eval_simpson (a, b, n, &func);
 233 
        comp = eval_simpson (a, b, n, &func);
217 
        abs = fabs(real - comp);
 234 
        abs = fabs(real - comp);
218 
        rel = fabs(real - comp) / fabs(real);
 235 
        rel = fabs(real - comp) / fabs(real);
219 
 
 236 
 
220 
        printf ("%.2d\t%e\t%e\t%e\t%e\n", n, h, comp, abs, rel);
 237 
        printf ("%.2d\t%e\t%e\t%e\t%e\n", i, h, comp, abs, rel);
221 
    }
 238 
    }
222 
 
 239 
 
223 
    return 0;
 240 
    return 0;
224 
}
 241 
}
225 
 
 242