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

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