MATLAB, second order differential equation, ode45

2차 미분방정식을 MATLAB 으로 바꾸어푸는 문제.

2차 미분방정식을 2개의 1차 미분방정식으로 변환시켰다! 

MATLAB 의 ode45 함수는 1차 미분방정식만 풀 수 있다.

        function dudt = odefun(t,u)

            dudt = [0;0];

            dudt(1) = u(2); % u1' = u2

            dudt(2) = -u(1)^n; % u2' = -u1^n

        end

odefun 에 두개의 1차 미분방정식을 넣고 코딩을 한다.

Solution:

clear; close;
figure(1)
for n = [1 5 9 15]
    [T,Y] = draw_oscillator(n);
    plot(T,Y)
    hold on
end
legend({'n = 1','n = 5','n = 9','n = 14'},'NumColumns',2)
xlabel('t')
ylabel('u(t)')
title("Solution of u'' + u^n = 0 ")
saveas(gcf,'oscillator','jpg')
figure(2)
for n = [1 5 9 15]
    [T,Y,Z] = draw_energy(n);
    Q = Z.^2 .* 0.5 + Y.^(n+1)/(n+1);  
    plot(T,Q)
    hold on
end
legend({'n = 1','n = 5','n = 9','n = 14'},'NumColumns',2)
xlabel('t')
ylabel('Kinetic + Potential Energy')
title("Total energy of u'' + u^n")
saveas(gcf,'energy','jpg')
function [T,Y] = draw_oscillator(n)
    u0 = [0 1];
    tspan = linspace(0,20);
        function dudt = odefun(t,u)
            dudt = [0;0];
            dudt(1) = u(2); % u1' = u2
            dudt(2) = -u(1)^n; % u2' = -u1^n
        end
    [T,U] = ode45(@odefun,tspan,u0);
    Y = U(:,1);
    end
function [T,Y,Z] = draw_energy(n)
    u0 = [0 1];
    tspan = linspace(0,20);
        function dudt = odefun(t,u)
            dudt = [0;0];
            dudt(1) = u(2); % u1' = u2
            dudt(2) = -u(1)^n; % u2' = -u1^n
        end
    [T,U] = ode45(@odefun,tspan,u0);
    Y = U(:,1);
    Z = U(:,2);
    end

Result:

Figure 2 에서 문제에 서술 된 것 처럼, n 이 작아질수록 0.5 에 가까워진다.

주어진 2개의 2차 미분방정식을 통해 projectile의 trail을 trace 하는 함수를 작성하자!

이번엔 ode45 말고 Euler's approximation method 를 이용해야한다.

Solution:

clear; close;
figure(1)
for n = [1 5 9 15]
    [T,Y] = draw_oscillator(n);
    plot(T,Y)
    hold on
end
legend({'n = 1','n = 5','n = 9','n = 14'},'NumColumns',2)
xlabel('t')
ylabel('u(t)')
title("Solution of u'' + u^n = 0 ")
saveas(gcf,'oscillator','jpg')
figure(2)
for n = [1 5 9 15]
    [T,Y,Z] = draw_energy(n);
    Q = Z.^2 .* 0.5 + Y.^(n+1)/(n+1);  
    plot(T,Q)
    hold on
end
legend({'n = 1','n = 5','n = 9','n = 14'},'NumColumns',2)
xlabel('t')
ylabel('Kinetic + Potential Energy')
title("Total energy of u'' + u^n")
saveas(gcf,'energy','jpg')
function [T,Y] = draw_oscillator(n)
    u0 = [0 1];
    tspan = linspace(0,20);
        function dudt = odefun(t,u)
            dudt = [0;0];
            dudt(1) = u(2); % u1' = u2
            dudt(2) = -u(1)^n; % u2' = -u1^n
        end
    [T,U] = ode45(@odefun,tspan,u0);
    Y = U(:,1);
    end
function [T,Y,Z] = draw_energy(n)
    u0 = [0 1];
    tspan = linspace(0,20);
        function dudt = odefun(t,u)
            dudt = [0;0];
            dudt(1) = u(2); % u1' = u2
            dudt(2) = -u(1)^n; % u2' = -u1^n
        end
    [T,U] = ode45(@odefun,tspan,u0);
    Y = U(:,1);
    Z = U(:,2);
    end

odeEuler 는 식을 받아 벡터로 결과를 저장한다.

Euler는 벡터를 받아 벡터를 결과로 저장한다.

X1에는 x''를 odeEuler로 풀어 x'의 값을 tspan에 따라 벡터형태로 저장되어있다.

X 에는 x'를 Euler로 풀어 x 의 값을 tspan에 따라 벡터형태로 저장되어있다.

y도 마찬가지.

Results:

>> [T,X,Y] = projectile(1,2,1,[0 5.9 0 0.8],0.01);

>> plot(X,Y)

>> title("Trace of projectile")

>> xlabel('x')

>> ylabel('y')

>> saveas(gcf,'trace','jpg')

(a)

Midpoint approximation을 MATLAB으로 구현한다.

Solution:

function

function [T,Y] = odeMidpoint(odefun,tspan,y0,h)
    T = tspan(1):h:tspan(2);
    if T(end) < tspan(2)
        T(length(T) + 1) = tspan(2);
    end
    Y = zeros(length(T),1);
    Y(1) = y0;
    for n = 1:length(T) - 2
        k1 = odefun(T(n),Y(n)); 
        k2 = odefun(T(n) + h/2 ,Y(n) + k1 * h/2); 
        Y(n+1) = Y(n) + k2 * h;    
    end
    
    h1 = T(end) - T(end - 1);
    k1 = odefun(T(end - 1),Y(end - 1)); 
    k2 = odefun(T(end - 1) + h1/2 ,Y(end - 1) + k1 * h1/2); 
    Y(end) = Y(end - 1) + k2 * h1;
end

(b) (c)

따라하면 된다..

clear; close;
odefun = @(t,y) -y^2 * t; y0 = 1;
fun = @(t) 2 / (t^2 + 2);
h = [1 0.5 0.2 0.1 0.01 0.001 0.0001];
tspan = [0 10];
error = zeros(length(h),1);

for i = 1:length(h)
    [T,Y] = odeMidpoint(odefun,tspan,y0,h(i));
    yN = Y(end);
    error(i) = abs(fun(10) - yN); % E(h)
    scatter(log(h(i)), log(error(i)));
    hold on;
end
grid on;
xlabel('log h'); ylabel('log E(h)');
legend('h=1','h=0.5','h=0.2','h=0.1','h=0.01','h=0.001','h=0.0001','Location','southeast','NumColumns',2)
title({"Global Truncation Error: Midpoint Method";"y' = -y^2t, y(0) = 1"})
saveas(gcf,'midpoint1','jpg');

clear; close;
odefun = @(t,y) -y^2 * t; y0 = 1;
fun = @(t) 2 / (t^2 + 2);
h = [1 0.5 0.2 0.1 0.01 0.001 0.0001];
tspan = [0 10];
error = zeros(length(h),1);

for i = 1:length(h)
    [T,Y] = odeMidpoint(odefun,tspan,y0,h(i));
    yN = Y(end);
    error(i) = abs(fun(10) - yN); % E(h)
    scatter(log(h(i)), log(error(i)));
    hold on;
end
grid on;
xlabel('log h'); ylabel('log E(h)');
legend('h=1','h=0.5','h=0.2','h=0.1','h=0.01','h=0.001','h=0.0001','Location','southeast','NumColumns',2)
title({"Global Truncation Error: Midpoint Method";"y' = -y^2t, y(0) = 1"})
saveas(gcf,'midpoint1','jpg');

Results:

(d)

polyfit 을 이용해 위 그래프에 맞는 slope를 찾는다.

Mathworks에 아주 설명이 잘 되어있다. https://www.mathworks.com/help/matlab/ref/polyfit.html

Solution:

p = polyfit(log(h),log(error'),1)

log - log 그래프에서 linear 함수의 slope 를 찾으므로 n = 1 이고, 이 결과값이 order 가 된다.