Algorithmization and Programming Examples from the 2-nd lecture ************************************************ Fixed point iterations (simple iterations) for solving x = f(x) ************************************************ Example 1.1 ============ f(x) = 0.2 * sin(x+0.5) Solving in the 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.m ------------- % script for simple iterations % - given the number of iterations % % f(x) = 0.2 * sin(x+0.5) x = input('give the starting point: '); n = input('give the number of iterations: '); for i=1:n x = 0.2 * sin(x+0.5); end disp('the result is'); disp(x); ------------- Ex1_1b.m ------------- % script for simple iterations, given the precision % % f(x) = 0.2 * sin(x+0.5) x_old = input('give the starting point: '); p = input('give the precision: '); x_new = 0.2 * sin(x_old + 0.5); while abs(x_old - x_new) > p x_old = x_new; x_new = 0.2 * sin(x_old + 0.5); end disp('the result is'); disp(x_new); ------------- Ex1_1c.m ------------- % script for simple iterations, given the precision, % with prevention of infinite cycle % % f(x) = 0.2 * sin(x+0.5) x_old = input('give the starting point: '); p = input('give the precision: '); c_max = input('give the maximal number of iterations: '); x_new = 0.2 * sin(x_old + 0.5); c = 1; % count of iterations while (abs(x_old - x_new) > p) & (c <= c_max) x_old = x_new; x_new = 0.2 * sin(x_old + 0.5); c = c + 1; end if c > c_max disp('number of iterations exceeded, the previous iteration is'); disp(x_old); end disp('the result is'); disp(x_new); ------------------------------ 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 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 >> Example 1.3 - the Jacobi method ================================= for solving a linear system A X = B ------------- Ex1_3.m ------------- % the Jacobi method - script for the conversion % of linear system A X = B to a system X = J X + V % D = diag(diag(A)); % diagonal of A L = tril(A, -1); % lower triangle of A (without the diagonal) P = triu(A, 1); % upper triangle of A (without the diagonal) U = - inv(D) * (L + P); % Jacobi iteration matrix V = inv(D) * B; % changed right hand side sp = max(abs(eig(U))); % spectral radius of J - criterion of convergence ---------------------------------- Solving in the Dialog window: >> A = [ 2 -1 0 ; -1 3 1 ; 0 1 2 ]; % system matrix >> B = [ 1 ; 3 ; 3 ]; % right hand side >> Ex1_3 % preparation for the Jacobi method >> X = V; % choose an initial value >> X = U * X + V % repeat, until X stays unchanged ... >> B - A * X % residual - should be close to zero ************************************************ Numerical integration ************************************************ Example 2.1 =========== f(x) = sin(sqrt(x)) , left Riemann sum ------------- Ex2_1.m ------------- % script for numerical integration of f(x) = sin(sqrt(x)) % on a given interval using left Riemann sum % a = input('give the beginning of the interval: '); b = input('give the end of the interval: '); n = input('give the number of integration points inside the interval: '); h = (b-a)/(n+1); % the integration step x = a:h:b; % the integration points fx = sin(sqrt(x)); % function values at the integration points S = h * sum( fx(1:n+1)); % left Riemann sum disp('the result is'); disp(S); ************************************************ 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 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 ****************************************************