Programming/MATLAB

MATLAB, implicit Euler Method, Global Truncation Error, BVPs, Fourier Transform

_Sudal 2019. 2. 16. 07:25

Implicit Euler Method 를 구현한다


function [T,Y] = odeIEuler(odefun,t,y0)
Y = zeros(length(t),1);
Y(1) = y0;
T = t;
    for n = 1:length(t) - 1
        h = T(n+1) - T(n);
        Y(n+1) = fzero(@(y) y - Y(n) - odefun(T(n+1),y) * h,0);
    end
end



 Forward, Backward and central difference formulas 를 이용해 Local Truncation Error 를 구현하는 script를 구현

%% (a)
clear; close;
fun = @(x) sin(sqrt(x));

x = 2;
k = 0:1:5;
fw = zeros(length(k),1);
bw = zeros(length(k),1);
c = zeros(length(k),1);

for k = k
   h = 10^(-k); 
   fw(k+1) = (fun(x + h) - fun(x)) / h; % forward method
   bw(k+1) = (fun(x) - fun(x - h)) / h; % backward method
   c(k+1) = (fun(x + h) - fun(x - h)) / (2*h); % central method
end


%% (b)
clear; close;
syms x
f(x) = sin(sqrt(x));
df(x) = diff(f(x),x);

k = 0:1:5;
fw = zeros(length(k),1);
bw = zeros(length(k),1);
c = zeros(length(k),1);
fwEh = zeros(length(k),1);
bwEh = zeros(length(k),1);
cEh = zeros(length(k),1);

h = 10.^(-k);
for k = k
    fw(k+1) = (double(f(2 + h(k+1))) - double(f(2))) / h(k+1); % forward method
    bw(k+1) = (double(f(2)) - double(f(2 - h(k+1)))) / h(k+1); % backward method
    c(k+1) = (double(f(2 + h(k+1))) - double(f(2 - h(k+1)))) / (2*h(k+1)); % central method

    fwEh(k+1) = abs(double(df(2)) - fw(k+1));
    bwEh(k+1) = abs(double(df(2)) - bw(k+1));
    cEh(k+1) = abs(double(df(2)) - c(k+1));
end

%% (c)
clear; close;
syms x
f(x) = sin(sqrt(x));
df(x) = diff(f(x),x);

k = 0:1:5;
fw = zeros(length(k),1);
bw = zeros(length(k),1);
c = zeros(length(k),1);
fwEh = zeros(length(k),1);
bwEh = zeros(length(k),1);
cEh = zeros(length(k),1);

h = 10.^(-k);
for k = k
    fw(k+1) = (double(f(2 + h(k+1))) - double(f(2))) / h(k+1); % forward method
    bw(k+1) = (double(f(2)) - double(f(2 - h(k+1)))) / h(k+1); % backward method
    c(k+1) = (double(f(2 + h(k+1))) - double(f(2 - h(k+1)))) / (2*h(k+1)); % central method

    fwEh(k+1) = abs(double(df(2)) - fw(k+1));
    bwEh(k+1) = abs(double(df(2)) - bw(k+1));
    cEh(k+1) = abs(double(df(2)) - c(k+1));
    
    subplot(3,1,1)
    scatter(log10(h(k+1)),log10(fwEh(k+1)));
    title('Absolute error for the forward difference formula for the first derivative')
    xlabel('log_1_0h')
    ylabel('log_1_0E(h)')
    ylim([-15 0])
    grid on    
    hold on
    
    subplot(3,1,2)
    scatter(log10(h(k+1)),log10(bwEh(k+1)));
    title('Absolute error for the backward difference formula for the first derivative')
    xlabel('log_1_0h')
    ylabel('log_1_0E(h)')
    ylim([-15 0])
    grid on
    hold on
    
    subplot(3,1,3)
    scatter(log10(h(k+1)),log10(cEh(k+1)));
    title('Absolute error for the central difference formula for the first derivative')
    xlabel('log_1_0h')
    ylabel('log_1_0E(h)')
    ylim([-15 0])

    grid on
    hold on
end

saveas(gcf,'fderrors1','jpg')

%% (d)
clear; close;
syms x
f(x) = sin(sqrt(x));
df(x) = diff(f(x),x);
d2f(x) = diff(df(x),x);

k = 0:1:5;
c = zeros(length(k),1);
cEh = zeros(length(k),1);

h = 10.^(-k);
for k = k
    c(k+1) = (double(f(2+h(k+1))) - 2*double(f(2)) + double(f(2-h(k+1))))/(h(k+1))^2;% central method    
    cEh(k+1) = abs(double(d2f(2)) - c(k+1));
    plot(log10(h(k+1)),log(cEh(k+1)),'o')
    hold on
    grid on
end
title('Absolute error for the central difference formula for the second derivative')
xlabel('log_1_0h')
ylabel('log_1_0E(h)')

saveas(gcf,'fderrors2','jpg')





second order BVPs 를 푸는 함수 구현


function [A,b] = bvpfd(p,q,r,xspan,BC,N)
    h = (xspan(2) - xspan(1))/(N+1);
    Q = ones(1,N); P1 = ones(1,N-1); P2 = ones(1,N-1);    
    b = ones(1,N);    
    for i = 1:1:N
        Q(i) = Q(i) * (2+q(i)*h^2);
    end
    
    b(1) = (p(1)*h/2+1)*BC(1) + (-p(1)*h/2+1)*BC(2) - h^2*r(1);
    if N > 1
        for i = 1:1:N-1
            P1(i) = P1(i) * (-p(i+1)*h/2 - 1);
            P2(i) = P2(i) * (p(i)*h/2 - 1);
        end
        
        b(1) = (p(1)*h/2 + 1)*BC(1) - h^2*r(1);
        b(end) = (-p(N)*h/2+1)*BC(2) - h^2*r(N);
        for i = 2:1:N-1
            b(i) = b(i) * (-h^2 * r(i));
        end
    end
    
    A = eye(N) .* Q + diag(P2,1) + diag(P1,-1);
    b = b';
end


ones(m,n): m * n 크기의 1 행렬을 만든다.

eye(N): N * N 크기의 행렬을 만든다. 대각선방향의 성분은 1로 채워진다

diag(M,n): M Matrix 를 n만큼 이동시켜 diagonal 하게 만들어준다. 나머지 항에는 0이 채워진다.


예)

>> A = [1 2 3 4 5];

>> diag(A,2)


ans =


     0     0     1     0     0     0     0

     0     0     0     2     0     0     0

     0     0     0     0     3     0     0

     0     0     0     0     0     4     0

     0     0     0     0     0     0     5

     0     0     0     0     0     0     0

     0     0     0     0     0     0     0


푸리에 급수를 반환하는 함수를 구현


function [X,Y] = partial_fourier(a0,a,b,L,x)

    X = x;

    Y = ones(1,length(X)) .* a0/2;

    for n = 1:1:length(a)

        Y = Y + a(n) .* cos(n*pi.*x/L) + b(n) .* sin(n*pi.*x/L);

    end

end


1 - Dimensinal Heat Equation에 관한 PDE를 푸는 함수 구현

function U = heat1D(k,b,L,x,t)
    p = length(t); q = length(x); N = length(b);
    U = zeros(p,q);
    for n = 1:1:N
        U = U + b(n) .* sin(n*pi.*x/L) .* exp(-n^2 * pi^2 / L^2 * k .* t');
    end
end