Image:Numerical integration illustration, h=1.png

I, the copyright holder of this work, hereby release it into the public domain. This applies worldwide. In case this is not legally possible: I grant anyone the right to use this work for any purpose, without any conditions, unless such conditions are required by law. Afrikaans | Alemannisch | Aragonés | العربية | Asturianu | Български | Català | Česky | Cymraeg | Dansk | Deutsch | Eʋegbe | Ελληνικά | English | Español | Esperanto | Euskara | Estremeñu | فارسی | Français | Galego | 한국어 | हिन्दी | Hrvatski | Ido | Bahasa Indonesia | Íslenska | Italiano | עברית | Kurdî / كوردی | Latina | Lietuvių | Latviešu | Magyar | Македонски | Bahasa Melayu | Nederlands | ‪Norsk (bokmål)‬ | ‪Norsk (nynorsk)‬ | 日本語 | Polski | Português | Ripoarisch | Română | Русский | Shqip | Slovenčina | Slovenščina | Српски / Srpski | Svenska | ไทย | Tagalog | Türkçe | Українська | Tiếng Việt | Walon | ‪中文(简体)‬ | ‪中文(繁體)‬ | zh-yue-hant | +/-

[edit] Matlab source code

% illustration of numerical integration
% compare the Forward Euler method, which is globally O(h) 
% with Midpoint method, which is globally O(h^2)
% and the exact solution

function main()

   f = inline ('y', 't', 'y'); % will solve y' = f(t, y)

   a=0; b=4; % endpoints of the interval where we will solve the ODE
   N = 5; T = linspace(a, b, N); h = T(2)-T(1); % the grid
   y0 = 1; % initial condition

   % solve the ODE
   Y_euler = solve_ODE (N, f, y0,  h, T, 1); % Forward Euler method
   Y_midpt = solve_ODE (N, f, y0,  h, T, 2); % midpoint method
   T_highres = a:0.1:b; Y_exact = exp(T_highres);
   
%  prepare the plotting window
   lw = 3; % curves linewidth
   fs = 20; % font size
   figure(1); clf; set(gca, 'fontsize', fs);   hold on;

   % colors
   red=[0.867 0.06 0.14];
   blue = [0, 129, 205]/256;
   green = [0, 200,  70]/256;

% plot the solutions
   plot(T, Y_euler, 'color', blue,  'linewidth', lw)
   plot(T, Y_midpt, 'color', green, 'linewidth', lw)
   plot(T_highres, Y_exact, 'color', red,   'linewidth', lw)

   % axes aspect ratio
   pbaspect([1 1.5 1]);

% save to disk
   disp(sprintf('Grid size is %0.9g', h))
   saveas(gcf, sprintf('Numerical_integration_illustration,_h=%0.2g.eps', h), 'psc2');
   
function Y = solve_ODE (N, f, y0,  h, T, method)

   Y = 0*T;
   
   Y(1)=y0;
   for i=1:(N-1)
          t = T(i); y = Y(i);

          if method == 1 % forward Euler method
                 
                 Y(i+1) = y + h*f(t, y);
                 
          elseif method == 2 % explicit one step midpoint method
                 
                 K = y + 0.5*h*f(t, y);
                 Y(i+1) =  y + h*f(t+h/2, K);
                 
          else
                 disp ('Don`t know this type of method');
                 return;
                 
          end
   end

• Numerical ordinary differential equations
• Midpoint method
• Euler method