# 數值分析 期末
[toc]
## Numerical differentiation
You can modified from the following codes into any other metohds and order of derivative.

%% Preparation
Y  = f(X);
dx = X(2)-X(1);

%% Second derivative with O(Δx^2)
function Yd = Diff1Ct(Y,dx)
    Yd          = zeros(1,length(X));
    Yd(1)       = (2*Y(1)-5*Y(2)+4*Y(3)-Y(4))/(dx^2);
    Yd(2:end-1) = (Y(1:end-2)-2*Y(2:end-1)+Y(3:end))/(dx^2);
    Yd(end)     = (2*Y(end)-5*Y(end-1)+4*Y(end-2)-Y(end-3))/(dx^2);
end

## Numerical integration
* Simpson's 1/3 method

    %% Preparation 
    Y = f(X);
    
    %% Simpson's 1/3 method
    function F = SolveIntS13(X,Y)
        dx = X(2)-X(1);
        F  = dx/3*(Y(1)+4*sum(Y(2:2:end))+2*sum(Y(3:2:end-1))+Y(end));
    end

* Simpson's 1/3 method

    %% Preparation 
    Y = f(X);
    
    %% trapezoidal method
    function F = SolveIntTrp(X,Y)
        dx = X(2)-X(1);
        F = dx/2*(Y(1)+Y(end))+dx*sum(Y(2:N)); 
    %   F = 0;
    %   for r = 1:length(X)-1
    %       F = F+(Y(r)+Y(r+1))*(X(r+1)-X(r))/2;
    %   end
    end

The commented codes (line 8-11) are for unevenly-spaced data.

## Solving ODE
### ODE
* Euler's method

    %% Preparation
    X = Xa:dx:Xb;
    
    %% Euler's Method
    function Y = SolveOdeEul(X,Y0,f)
        dx   = X(2)-X(1);
        Y    = zeros(size(X));
        Y(1) = Y0;
        for i = 1:length(X)-1
            k      = f(X(i),Y(i));
            Y(i+1) = Y(i)+k*dx;
        end
    end

* Modified Euler's method

    %% Preparation
    X = Xa:dx:Xb;
    
    %% Modified Euler's Method
    function Y = SolveOdeEuM(X,Y0,f)
        dx   = X(2)-X(1);
        Y    = zeros(size(X));
        Y(1) = Y0;
        for i = 1:length(X)-1
            k1     = f(X(i),Y(i));
            Y(i+1) = Y(i)+k1*dx;
            k2     = f(X(i+1),Y(i+1));
            Y(i+1) = Y(i)+(k1+k2)/2*dx;
        end
    end

* 4th Order Runge-Kutta Method (RK4)

    %% Preparation
    X = Xa:dx:Xb;
    
    %% 4th Order Runge-Kutta Method
    function Y = SolveOdeRK4(X,Y0,f)
        dx   = X(2)-X(1);
        Y    = zeros(size(X));
        Y(1) = Y0;
    for i = 1:length(X)-1
            k1     = f(X(i),Y(i));
            k2     = f(X(i)+dx/2,Y(i)+k1*dx/2);
            k3     = f(X(i)+dx/2,Y(i)+k2*dx/2);
            k4     = f(X(i)+dx,Y(i)+k3*dx);
            Y(i+1) = Y(i)+(k1+2*k2+2*k3+k4)/6*dx;
        end
    end

### ODE system
* Euler's method

    %% Preparation
    X = Xa:dx:Xb;
    
    %% Euler's Method
    function [Y1,Y2] = SolveOdeSetEul(X,f1,Y01,f2,Y02)
        dx = X(2)-X(1);
        Y1 = zeros(size(X)); Y1(1) = Y01;
        Y2 = zeros(size(X)); Y2(1) = Y02;
        k  = [0,0];
        for i = 1:length(X)-1
            k(1)  = f1(Y1(i),Y2(i));
            k(2)   = f2(Y1(i),Y2(i));
            Y1(i+1) = Y1(i)+k(1)*dx;
            Y2(i+1) = Y2(i)+k(2)*dx;
        end
    end

* 4th Order Runge-Kutta Method

    %% Preparation
    X = Xa:dx:Xb;
    
    %% 4th Order Runge-Kutta Method
    function [Y1,Y2] = SolveOdeSetRK4(X,f1,Y01,f2,Y02)
        dx = X(2)-X(1);
        Y1 = zeros(size(X)); Y1(1) = Y01;
        Y2 = zeros(size(X)); Y2(1) = Y02;
        for i = 1:length(X)-1
            k1(1)   = f1(Y1(i),Y2(i));
            k1(2)   = f2(Y1(i),Y2(i));
            k2(1)   = f1(Y1(i)+k1(1)*dx/2,Y2(i)+k1(2)*dx/2);
            k2(2)   = f2(Y1(i)+k1(1)*dx/2,Y2(i)+k1(2)*dx/2);
            k3(1)   = f1(Y1(i)+k2(1)*dx/2,Y2(i)+k2(2)*dx/2);
            k3(2)   = f2(Y1(i)+k2(1)*dx/2,Y2(i)+k2(2)*dx/2);
            k4(1)   = f1(Y1(i)+k3(1)*dx,Y2(i)+k3(2)*dx);
            k4(2)   = f2(Y1(i)+k3(1)*dx,Y2(i)+k3(2)*dx);
            Y1(i+1) = Y1(i)+(k1(1)+2*k2(1)+2*k3(1)+k4(1))/6*dx;
            Y2(i+1) = Y2(i)+(k1(2)+2*k2(2)+2*k3(2)+k4(2))/6*dx;
        end
    end

## 考古
1. Consider a stratified lake with a surface area $A=4 \,\mathrm{km}^2$ with two well mixed layers of thickness $h_1=1 \,\mathrm{m}$ and $h_2=10 \,\mathrm{m}$. A river flows into the top layer at a rate $Q=12\cdot10^6 \,\mathrm{m}^3\cdot\mathrm{yr}^{-1}$ and another river flows out of the top layer at an equal rate $Q$.
* Suppose that a gold mine releases some element into the river flowing into the lake such that the concentration in the river is $C_i=1 \,\mathrm{mol}\cdot\mathrm{m}^{-3}$.
* Suppose furthermore that turbulent mixing in the lake causes an exchange of volume $q=12\cdot10^4 \,\mathrm{m}^3\cdot\mathrm{yr}^{-1}$ between the layers.
* Assuming there are no other sources or sinks of the element.
1. Set up a system of differential equations describing the time evolution of the concentrations ($C_1$, $C_2$) of the tracer in the upper and lower layers of the lake.
\begin{align}
\frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix}
C_1 \\ C_2\end{bmatrix} =
\frac{1}{A} \left(
\begin{bmatrix}
-(Q+q)/h_1 & q/h_1 \\ q/h_2 & -q/h_2
\end{bmatrix}\begin{bmatrix}
C_1 \\ C_2\end{bmatrix} +\frac{Q}{h_1}
\begin{bmatrix}C_i \\ 0\end{bmatrix}
\right).
\end{align}
2. Compute the equilibrium concentrations $C_1$ and $C_2$ at steady state.
\begin{align}C_1=C_2=C_i.\end{align}
4. Suppose that the mine suddenly shuts down and stops polluting. How long will it take before $C_2$ drops to a level that is $1\%$ of its equilibrium concentration during the time that the mine was operating? Please use the Euler Forward method to numerically integrate your system of equations. Estimate the minimum number of years that would be needed. You do not need to worry about the accuracy of the estimate, but you must use a time step that is small enough for the Euler Forward integration method to be stable.

    dt      = 1/12;
    T       = 0:dt:10000;
    [C1,C2] = SolveOdeSetEul(T,@F1,1,@F2,1);
    [~,t]   = max(C2<=C2(1)*0.01);
    t       = T(t);

    function F = F1(C1,C2)
        A  = 4e6;
        h1 = 1;
        Q  = 12e6;
        q  = 12e4;
        Ci = 0;
        F  = 1/A*(-(Q+q)/h1*C1+q/h1*C2+Q/h1*Ci);
    end
    function F = F2(C1,C2)
        A  = 4e6;
        h2 = 10;
        q  = 12e4;
        F  = 1/A*(q/h2*C1-q/h2*C2);
    end

$t\sim1550\,\mathrm{yr}$.