nul, your issue is about understanding the mathematical method, not so much a programming issue.
That said, perhaps the following variant of the same code helps you undestand better:
Code:
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
/* Slope function f(t, x) = y'(t) where x = y(t) */
static double f(const double t, const double x)
{
const double tt = t * t;
return 2.0 * t * tt * (tt / (x * x)); /* 2 * t^5 / x^2 */
}
/* Known exact function y(t) */
static double known_y(const double t)
{
const double t2 = t * t;
const double t4 = t2 * t2;
return cbrt(t2 * t4 - 37.0);
}
int main(int argc, char *argv[])
{
const double t0 = 2.0; /* Starting point */
const double y0 = 3.0; /* Function at starting point */
const double t1 = 3.0; /* Ending point */
double h, y;
long steps, i;
char dummy;
if (argc != 2 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) {
fprintf(stderr, "\n");
fprintf(stderr, "Usage: %s [ -h | --help ]\n", argv[0]);
fprintf(stderr, " %s STEPS\n", argv[0]);
fprintf(stderr, "\n");
return EXIT_FAILURE;
}
if (sscanf(argv[1], " %ld %c", &steps, &dummy) != 1 || steps < 1L) {
fprintf(stderr, "%s: Invalid number of steps.\n", argv[1]);
return EXIT_FAILURE;
}
h = (t1 - t0) / (double)steps;
y = y0;
printf("# t f(t,y) y y(t) h\n%.9f %.9f %.9f %.9f %.9f\n", t0, f(t0, y0), y0, known_y(t0), h);
for (i = 1L; i <= steps; i++) {
const double phase = (double)i / (double)steps;
const double t = (1.0 - phase) * t0 + phase * t1;
const double fvalue = f(t, y);
y += h * fvalue;
printf("%.9f %.9f %.9f %.9f\n", t, fvalue, y, known_y(t));
}
fflush(stdout);
fprintf(stderr, "Expected y = %.16f\n", known_y(t1));
fprintf(stderr, "Evaluated y = %.16f\n", y);
return EXIT_SUCCESS;
}
The Wikipedia Euler method article has a pretty good description. Let me re-describe it, with variable names matching the above program and the article.
The problem at hand is an initial value problem. We know that
dy(t) / dt = f(t, y(t))
f(t, x) = 2 t5 / x2
y(2) = 3.0
and we want to find out
y(3) = ?
In other words, we don't know the function y(t), but we know its slope is f(t,x) where x=y(t). One way to approximate this is to apply the Euler method, where t=t0+nh, and h is the number of steps between t0=2 (the starting point) and t1=3 (the ending point), and we know the initial value y(t0)=y0=3.
The above program has hardcoded the starting point (t0=2.0, t1=3.0, and y0=3.0), but the number of steps taken is supplied on the command line. The step size, h, is determined from those three (h = (t1 - t0) / steps).
Because we have an integer number of steps (iterations), I'm using a stable numeric trick to compute t at each step. The number of iterations (starting at 1) divided by the total number of steps gives us phase, a value beween 0 (exclusive) and 1 (inclusive), at very high precision. We can then calculate t = (1.0-phase)*t0 + phase*t1 to calculate t without cumulative errors. (Here, the difference to just adding h to t at each step would not make much of a difference, but if the start and end point were dozens of orders of magnitude different, then this method shines: it will yield reliable values even then.)
During each iteration, in the loop body, we precalculate t and fvalue = f(t, y), then update y=y+h*fvalue. This is equivalent to the Wikipedia article yN+1=yN+h×f(tN,yN). We use separate variables, so we can print the actual values used, as well as yN+1, the point of the iteration.
After the loop is done, we have the approximation of y(t1) in y. To make it easier to see the effect the number of steps has on the approximation, the summary is printed to standard error (while all the intermediate values are printed to standard output). Note that the loop will run much faster, if only the end result were to be printed; printing floating-point values is kind-of slow, when you have lots of iterations.
I compiled and ran the above as euler, with different numbers of steps. The results, summarized, are as follows:
Code:
# Steps Approximation h
1 57.0000000000000000 1
2 15.1173635453141522 0.5
3 11.2860983089115550 0.3333333
5 10.0388388035749934 0.2
10 9.4014859985922055 0.1
20 9.1159639917882256 0.05
50 8.9518255220655405 0.02
100 8.8981953539855372 0.01
500 8.8556663104274111 0.002
1000 8.8503733129637165 0.001
10000 8.8456139814468280 0.0001
100000 8.8451382755005721 0.00001
1000000 8.8450907071763734 0.000001
10000000 8.8450859503674462 0.0000001
# Exact 8.8450854218326640
nul, in your code delta_t = .01 = h, which translates to just 100 steps. As the Wikipedia article suggests, the error is roughly proportional to the step size, so you only get about two digits of precision. Using more steps (more values of t and smaller h (your delta_t) gives you more precision.
Questions?
Originally Posted by nul
