Algorithmization and Programming Examples from the 2-nd lecture rewritten to C ************************************************ Fixed point iterations (simple iterations) for solving x = f(x) ************************************************ Example 1.1 ============ f(x) = 0.2 * sin(x+0.5) Solving in the Matlab Dialog window: >> x=0.4; >> x=0.2*sin(x+0.5) x = 0.15667 >> x=0.2*sin(x+0.5) x = 0.12210 >> ... 0.11655 ... 0.11564 ... 0.11550 x = 0.11547 >> x=0.2*sin(x+0.5) x = 0.11547 >> ------------- Ex1_1a.c ------------- #include #include int main() { /* program for simple iterations - given the number of iterations f(x) = 0.2 * sin(x+0.5) */ int n, i; float x; printf("give the starting point: "); scanf("%f",&x); printf("give the number of iterations: "); scanf("%d",&n); for(i=1; i<=n; i++) { x = 0.2 * sin(x+0.5); } printf("the result is %f\n", x); return 0; } ------------- Ex1_1b.c ------------- #include #include int main() { /* program for simple iterations, given the precision f(x) = 0.2 * sin(x+0.5) */ float x_old, x_new, p; printf("give the starting point: "); scanf("%f",&x_old); printf("give the precision: "); scanf("%f",&p); x_new = 0.2 * sin(x_old + 0.5); while ( fabs(x_old - x_new) > p ) { x_old = x_new; x_new = 0.2 * sin(x_old + 0.5); } printf("the result is %f\n", x_new); return 0; } ------------- Ex1_1c.c ------------- #include #include int main() { /* program for simple iterations, given the precision with prevention of infinite cycle f(x) = 0.2 * sin(x+0.5) */ int c_max, c; float x_old, x_new, p; printf("give the starting point: "); scanf("%f",&x_old); printf("give the precision: "); scanf("%f",&p); printf("give the maximal number of iterations: "); scanf("%d",&c_max); x_new = 0.2 * sin(x_old + 0.5); c = 1; // count of iterations while ( (fabs(x_old - x_new) > p) && (c <= c_max) ) { x_old = x_new; x_new = 0.2 * sin(x_old + 0.5); c = c + 1; } if (c > c_max) { printf("number of iterations exceeded, the previous iteration is %f,\n", x_old); } printf("the result is %f\n", x_new); return 0; } ------------------------------ Example 1.2 - logistic map ========================== (discrete time demographic model - reproduction vs. starvation) f(x) = r*x*(1-x) where r ... parameter x ... population size ... in < 0, 1> a fixed point (steady state) exists for r in ( 1, 3) Solving in the Matlab Dialog window: >> r=2; >> x = 0.2; % any number in < 0, 1> >> x = r*x*(1-x) x = 0.32000 >> x = r*x*(1-x) x = 0.43520 >> ... 0.49160 ... 0.49986 x = 0.50000 >> x = r*x*(1-x) x = 0.50000 >> ------------- Ex1_2.c ------------- #include int main() { /* prints fixed-point iterations for a logistic map f(x) = r*x*(1-x) input: r ... parameter of the map x ... starting value of the population size n ... number of iterations */ int n, i; float r, x; printf("give the parameter, in ( 0, 3): "); scanf("%f",&r); printf("give the number of iterations: "); scanf("%d",&n); printf("give the starting point, in < 0, 1>: "); scanf("%f",&x); printf("%10f\n", x); for(i=1; i<=n; i++) { x = r*x*(1-x); printf("%10f\n", x); } return 0; } ************************************************ Numerical integration ************************************************ Example 2.1 =========== f(x) = sin(sqrt(x)) , left Riemann sum ------------- Ex2_1.c ------------- #include #include int main() { /* script for numerical integration of f(x) = sin(sqrt(x)) on a given interval using left Riemann sum */ int n, k; float a, b, h, x, fx, S; printf("give the beginning of the interval: "); scanf("%f",&a); printf("give the end of the interval: "); scanf("%f",&b); printf("give the number of integration points inside the interval: "); scanf("%d",&n); h = (b-a)/(n+1); // the integration step S = 0; // left Riemann sum for (k=0; k <= n; k++) { x = a + k*h; // integration point fx = sin(sqrt(x)); // function value at the integration point S = S + fx; } S = h * S; printf("the result is %f\n", S); return 0; } ************************************************ Explicit Euler's method for solving y' = f(x, y) ************************************************ Example 3.1 ============ f(x) = y / x^2, x_0 = 1, y_0 = 2 Solving in the Matlab Dialog window: >> h = 0.2; % choose a step-size >> x = 1; y = 2; % initial condition - starting point Repeat the following 4 steps: >> f = y / x^2; % f ... direction of the tangent line at [x, y] >> x = x + h; % new value of x >> y = y + f * h; % new value of y (on the tangent line) >> disp(y); % print the new value on the screen ------------- Ex3_1.c ------------- #include #include int main() { /* Explicit Euler's method for solving y' = y / x^2 */ int n, i; float h, x, y, f; printf("give the initial value of x: "); scanf("%f",&x); printf("give the initial value of y: "); scanf("%f",&y); printf("give the step-size: "); scanf("%f",&h); printf("give the number of steps: "); scanf("%d",&n); printf("%10.2f %10.4f\n", x, y); // print the starting point for(i=1; i<=n; i++) { f = y / pow(x,2); // f ... direction of the tangent line at [x, y] x = x + h; // new value of x y = y + f * h; // new value of y (on the tangent line) printf("%10.2f %10.4f\n", x, y); // print the new value of [x, y] } return 0; } ****************************************************