Algorithmization and Programming

=======================================
=========    FUNCTIONS    =============
=======================================

 ... advanced


=====================================================
=====================================================
Passing a function as a parameter to another function
=====================================================

### Example 5

************************************************
Explicit Euler's method
for solving  y' = f(x, y)
************************************************


  f(x,y) = y / x^2,  x_0 = 1,  y_0 = 2


===================
==  file Eul.m  ===
===================

function y = Eul(f, x0, y0, h, n)
%
% input: 
%   f(x,y)  - link to the rhs of the differential eq.
%   x0, y0  - initial condition (starting point)
%   h       - a step-size
%   n       - number of steps 

  x = x0;
  y = zeros(1,n);
  y(1) = y0;
  for k=2:n
    y(k) = y(k-1) + f(x,y(k-1)) * h; 
    x = x + h;       
  end
end


=================
==  file f.m  ===
=================

function z = f(x,y)
% f(x,y) = y / x^2; 
  z = y / x^2; 
end

---------------

  --- calling the function Eul:

>> p = Eul(@f,1,2,0.2,5)     % @f = Function Handle


================================================


### Example 6

function S = integ(f, a, b, n)
% 
% numerical integration of f on (a, b)
%
% input: 
%   f ... link to a function of one parameter
%   a , b ... endpoints of the interval
%   n ... the number of subintervals 
%
% output:
%  S ... integral from a to b of f
%
  x = linspace(a,b,n+1); % the integration points
  h = x(2)-x(1);         % the integration step
  fx = f(x);             % function values at integration points
  S = h * sum( fx(1:n)); % left Riemann sum
end

==================

>> integ( @sin, 0, pi, 20)  % @sin = Function Handle

ans =

    1.9959

>> g = @(x) x.^2 - 3*x      % g = Function Handle

g = 

    @(x)x.^2-3*x

>> integ(g, 1, 3, 20 )  % exact value: -3.333...

ans =

   -3.4300

>>

==============================================


GLOBAL VARIABLES

>> global g = 9.81;
>> t = fall(10)

====================
==  file fall.m  ===
====================

function t = fall(h)
% computes time it takes for the object to drop
% input: h - height of a tower [m]

  global g   % gravitational acceleration
  t = sqrt(2*h/g);
end


===================
=== some script ===
===================

global g = 9.81;
global r = 1000;

p = press(10);


=====================
==  file press.m  ===
=====================

function p = press(h)
% computes hydrostatic pressure p = h*g*r  [m]*[N/kg]*[kg/m3]
% input: h - height of water column [m]

  global g   % gravitational acceleration
  global r   % the specific weight of water

  p  = h*g*r;
end

====================================
====================================