
% funktionsnamnen nedan h?nvisar till tex Edsberg upg 7.1 (ed7p1) eller
% Timothy Sauer (sauer5p11)

function ovn121115

clc, clear all, close all

%ed6p1 % Ta bort kommentar för att köra uppgiften
ed7p1

function sauer5p11

% Hitta en differensformel f?r f'(x) i punkterna f(x-h), f(x), f(x+3h)
% L?sning: det g?ller att v?lja a,b,c s?dana att termen f?r f(x) och f''(x) elimineras i
% taylorutvecklingen, och att f'(x) blir kvar + h?gre ord. feltermer.

% f'(x) = (a*f(x-h) + b*f(x) + c*f(x+3h)/h

% f(x+h) = f(x) + f'(x)*h + f''(x)/2*h^2 + f'''(x)/6*h^3

% f'(x) = a*(f(x) - f'(x)*h + f''(x)/2*h^2 - f'''/6*h^3) +
%         b*(f(x)) +
%         c*(f(x) + f'(x)*3h + f''(x)*9/2*h^2 - f'''*27/6*h^3) 

% Man f?r d? ekvationssystemet: 
% A*[a b c]' = [0 1 0]', d?r 0 representerar en elimnerad term, och etta
% en ej eliminerad term, och d?r 

A = [1 1 1; -1 0 3; 1 0 9]

% och l?sningen blir d?

A\[0 1 0]'

% och dominerande felterm blir O(h^3)/h = O(h^2)



function ed7p1
% L?sa en diff-ekvation med Euler fram?t och j?mf?r med exakta l?sningen.
% @ ul?ses "funktion av"

yprim = @(t, y) 1 + t - y; % Diff-ekvationen
yex = @(t) exp(-t) + t; % Exakta l?sningen

yex(0.8)

for h = [0.05 0.1 0.2] % F?r vart och ett av dess h-v?rden g?r det som st?r i slingan
   y = 1;
   
   for t = 0:h:(0.8-h) % Ber?kna med Euler fr?n t = 0 till 0.8. N?r man st?r
       % i t = 0.8-h skall man ju ta sista steget f?r att komma till t = 0.8, d?rav "(0.8)-h".
      y  = y + yprim(t, y)*h;
      t
   end
   y
end


function ed6p1
% Ber?kna n?gra integraler numeriskt, d?r en ?r osn?ll.

f = @(x) exp(x)./(1+2*x.^3);

% Ber?kna integral av f med en viss tolerans. Testa att ?ndra toleransen och se
% hur resultatet p?verkas.
integral(f, 0, 3, 'AbsTol', 1e-6, 'RelTol', 0.5e-4)

g = @(x) sin(x)./x; % Osn?ll integrand som 'integral' dock ändå klarav av.
 
integral(@gg, 0, pi) % Integrera den snälla integranden.

fplot(f, [0 3])


% Monte-Carlo-integration av f ovan. F?r stora N skall det ge samma svar som
% 'integral' ovan.
N = 100000;

x = rand(1,N) * 3; % Slumpa punkter i intervallet.
% x = linspace(0, 3, N); % J?mf?r med vanlig Rieman-summa

ff = f(x);

sum(ff)*(3/N) % summera ihop alla punkter och multiplicera med "steget" (3/N)

function ggut = gg(x) % Sn?ll integrand.
if x == 0
    ggut = 1;
else 
    ggut = sin(x)./x;
end



