Lab 8 - Kamil Breguła

function [t, y] = euler(f, tspan, y0, N)
    yy = [y0];
    xx = [0];

    a = tspan(1);
    b = tspan(2);
    h = (b - a) / N;

    for i = 2:N+1
        current_x = xx(i-1) + h;
        current_y = yy(i-1) + h * f(current_x, yy(i-1));
        xx = [xx current_x];
        yy = [yy current_y];
        if current_x > 20
            break;
        end;
    end
    t = xx;
    y = yy;

function [t, y] = euler_m(f, tspan, y0, N)
    yy = [y0];
    xx = [0];

    a = tspan(1);
    b = tspan(2);
    h = (b-a)/N;
    for i = 2:N+1
        x_i = xx(i-1);
        x_i_1 = xx(i-1) + h;

        y_i = yy(i-1)
        
        y_i_1 = y_i + h * f(x_i + (1/2) * h, y_i + 1/2 * h * f(x_i, y_i))
    
        xx = [xx x_i_1];
        yy = [yy y_i_1];
        
        if x_i_1 > 20
            break;
        end;
    end
    t = xx;
    y = yy;

function [x,y] = ab3(f, tspan, y0, n)
    a = tspan(1);
    b = tspan(2);
    
    h = (b - a) / n;
    x = (a : h : b);
    
    y(1) = y0;
    
    for i = 1 : 2
        k_1 = h * feval(f, x(i), y(i));
        k_2 = h * feval(f, x(i) + h / 2, y(i) + k_1/2);
        
        y(i+1) = y(i) + k_2;
    end
    
    z0 = feval(f, a, y0);
    z(1) = feval(f, x(1), y(1));
    z(2) = feval(f, x(2), y(2));
    
    y(3) = y(2) + h / 12 * (23 * z(2) - 16 * z(1) + 5 * z0);
    
    for i = 3 : n
        z(i) = feval(f, x(i), y(i));
        y(i+1) = y(i) + h / 12 * (23 * z(i) - 16 * z(i-1) + 5 * z(i-2));
    end

clear;
clf;
cla;

syms y(x)
dsolve(diff(y)/diff(x) - x - y, y(0) == 1)
dsolve(x^2*diff(x)-y*diff(x)-diff(y),  y(0) == 1)
dsolve(diff(y) == x^2 - y, y(0) == 1)


fn1 = @(x, y) x+y
fn2 = @(x, y) x.^2 * x - y * x
fd1 = @(x) 2 * exp(x) - x - 1 
fd1 = @(x) 2 * exp(x) - x - 1
fd2 = @(x) x.^2 - exp(-x) - 2 * x + 2
fd2 = @(x) x.^2 - exp(-x) - 2 * x + 2
a = 0
b = 10
tspan = [a b];

precision = 0.01

[x1a, y1a] = euler(fn1, tspan, 1, (b - a) / precision)
[x1b, y1b] = euler(fn2, tspan, 1, (b - a) / precision)

[x2a, y2a] = euler_m(fn1, tspan, 1, (b - a) / precision)
[x2b, y2b] = euler_m(fn2, tspan, 1, (b - a) / precision)

[x3a, y3a] = ab3(fn1, tspan, 1, (b-a) / precision)
[x3b, y3b] = ab3(fn2, tspan, 1, (b-a) / precision)

x4a = [0:precision:10]
y4a = fd1(x4a)

x4b = [0:precision:10]
y4b = fd2(x4b)

figure(1)
subplot(3, 2, 1)
title('Wartosci funkcji precision = 0.01')
hold on
plot(x1a, y1a)
plot(x2a, y2a)
plot(x3a, y3a)
plot(x4a, y4a)
legend('euler', 'euler zmodyfikowawny', 'ab3', 'dsolve')
hold off

subplot(3, 2, 2)

err1 = (fd1(x1a) - y1a).^2
err2 = (fd1(x2a) - y2a).^2
err3 = (fd1(x3a) - y3a).^2
err1_p1_x = x1a
err2_p1_x = x2a
err3_p1_x = x3a
err1_p1_y = err1
err2_p1_y = err2
err3_p1_y = err3

title("Wartosc bledu precission = 0.01")
hold on
plot(x1a, err1)
plot(x2a, err2)
plot(x3a, err3)
legend('euler', 'euler zmodyfikowany', 'ab3')
hold off

precision = 0.1

[x1a, y1a] = euler(fn1, tspan, 1, (b - a) / precision)
[x1b, y1b] = euler(fn2, tspan, 1, (b - a) / precision)

[x2a, y2a] = euler_m(fn1, tspan, 1, (b - a) / precision)
[x2b, y2b] = euler_m(fn2, tspan, 1, (b - a) / precision)

[x3a, y3a] = ab3(fn1, tspan, 1, (b-a) / precision)
[x3b, y3b] = ab3(fn2, tspan, 1, (b-a) / precision)

x4a = [0:precision:10]
y4a = fd1(x4a)

x4b = [0:precision:10]
y4b = fd2(x4b)

subplot(3, 2, 3)
title('Wartosci funkcji precision = 0.1')
hold on
plot(x1a, y1a)
plot(x2a, y2a)
plot(x3a, y3a)
plot(x4a, y4a)
legend('euler', 'euler zmodyfikowawny', 'ab3', 'dsolve')
hold off

subplot(3, 2, 4)

err1 = (fd1(x1a) - y1a).^2
err2 = (fd1(x2a) - y2a).^2
err3 = (fd1(x3a) - y3a).^2
err1_p2_x = x1a
err2_p2_x = x2a
err3_p2_x = x3a
err1_p2_y = err1
err2_p2_y = err2
err3_p2_y = err3

title("Wartosc bledu precission = 0.1")
hold on
plot(x1a, err1)
plot(x2a, err2)
plot(x3a, err3)
legend('euler', 'euler zmodyfikowany', 'ab3')
hold off

precision = 0.5

[x1a, y1a] = euler(fn1, tspan, 1, (b - a) / precision)
[x1b, y1b] = euler(fn2, tspan, 1, (b - a) / precision)

[x2a, y2a] = euler_m(fn1, tspan, 1, (b - a) / precision)
[x2b, y2b] = euler_m(fn2, tspan, 1, (b - a) / precision)

[x3a, y3a] = ab3(fn1, tspan, 1, (b-a) / precision)
[x3b, y3b] = ab3(fn2, tspan, 1, (b-a) / precision)

x4a = [0:precision:10]
y4a = fd1(x4a)

x4b = [0:precision:10]
y4b = fd2(x4b)

subplot(3, 2, 5)
title('Wartosci funkcji precision = 0.5')
hold on
plot(x1a, y1a)
plot(x2a, y2a)
plot(x3a, y3a)
plot(x4a, y4a)
legend('euler', 'euler zmodyfikowawny', 'ab3', 'dsolve')
hold off

subplot(3, 2, 6)

err1 = (fd1(x1a) - y1a).^2
err2 = (fd1(x2a) - y2a).^2
err3 = (fd1(x3a) - y3a).^2
err1_p3_x = x1a
err2_p3_x = x2a
err3_p3_x = x3a
err1_p3_y = err1
err2_p3_y = err2
err3_p3_y = err3

title("Wartosc bledu precission = 0.5")
hold on
plot(x1a, err1)
plot(x2a, err2)
plot(x3a, err3)
legend('euler', 'euler zmodyfikowany', 'ab3')
hold off


figure(2)

title("Wartosc bledu precission = 0.01, 0.1, 0.5")
hold on
plot(err1_p1_x, err1_p1_y, 'x-')
plot(err2_p1_x, err2_p1_y, 'x-')
plot(err3_p1_x, err3_p1_y, 'x-')
plot(err1_p2_x, err1_p2_y, '+-')
plot(err2_p2_x, err2_p2_y, '+-')
plot(err3_p2_x, err3_p2_y, '+-')
plot(err1_p3_x, err1_p3_y, 'd-')
plot(err2_p3_x, err2_p3_y, 'd-')
plot(err3_p3_x, err3_p3_y, 'd-')

legend( ...
    'euler precision=0.01', 'euler zmodyfikowany precision=0.01', 'ab3 precision=0.01', ...
    'euler precision=0.1', 'euler zmodyfikowany precision=0.1', 'ab3 precision=0.1', ...
    'euler precision=0.5', 'euler zmodyfikowany precision=0.5', 'ab3 precision=0.5')
hold off


Widok normalny

Zbliżenie I

Zbliżenie II