---
tags: homeworks, Numerical Methods
---
# Numerical Methods Homework-8
## B10702026 林琨霖
[Online link of this homework](https://hackmd.io/4pJSURBoScmHiDbOTi6vpQ?view)
It's recommanded to review this homework through the link above. The content is same with pdf but the online link is with a better layout, clickable index and syntax highlighting.
> For the online link, you can click `...` at the upper right corner and then click `版本與Github同步` to check for the last edition time, making sure this homework is finished before the dead line.
> 
>
> 
---
## Q1 Evaluate the following integral:
$$
\int _0 ^{\pi /2}(8+4\text{ cos }x)dx
$$
For each of the numerical estimates (b) through (g), determine the true percent
relative error based on (a).
### (a) Analytically
$$
\begin{aligned}
\int _0 ^{\pi /2}(8+4\text{ cos }x)dx &=(8x+4\text{sin} x) \vert_0 ^{\pi /2} \\
&= 4\pi+4 \\
&= 16.566370614359172
\end{aligned}
$$
### (b) Single application of the trapezoidal rule
- Code

  f = @(x) 8 + 4.*cos(x); %// function to integral
  from = 0;               %// upper limit of integral
  to   = pi/2;            %// lower limit of integral
  ANS  = (8*to + 4*sin(to)) - (8*from + 4*sin(from));

  %// Single application of the Trapezoidal Rule
  I = (to-from) * (f(from) + f(to)) / 2;

  %// Print the result -----------------------------
  fprintf("Single application of the Trapezoidal Rule.\n");
  fprintf("Estiamted Integral:          %.15f\n", I);
  fprintf("True value:                  %.15f\n", ANS);
  fprintf("True relative percent error: %.15f%%\n\n\n", 100*abs( (ANS-I)/ANS ));

  %// Plot -----------------------------------------
  x = linspace(from-0.5, to+0.5, 100);
  hold on;
  plot(x, f(x));
  plot([from from to to], [0 f(from) f(to) 0], '-');
  legend("Original function", "Approximated Trapezoid");
  title("Single Trapezoidal Rule");

- Result

  >> HW08_Q1_b
  Single application of the Trapezoidal Rule.
  Estiamted Integral:          15.707963267948966
  True value:                  16.566370614359172
  True relative percent error: 5.181625875652982%

I plot the original function and the approximated one with trapezoidal rule.
The estimated integral is the area closed by the red lines.
The function to integral $8+4\text{ cos }x$ is a concave function and the red trapezoid obviously under-estimates the true integral.
\-

### (c\) Composite trapezoidal rule with $n=2$ and $n=4$
- Code

  %// Declartions
  f = @(x, varargin) 8 + 4.*cos(x); %// function to integral
  from = 0;                         %// upper limit of integral
  to   = pi/2;                      %// lower limit of integral
  n    = 4;                         %// how many segaments
  ANS  = (8*to + 4*sin(to)) - (8*from + 4*sin(from));


  %// Composite Trapezoidal Rule
  step = (to-from) / n; %// aka h, evenly spaced interval
  x = linspace(from+step, to-step, n-1);
  I = step/2 * (f(from) + f(to) + 2*sum(f(x)));

  %// Print the result -----------------------------
  fprintf("Composite Trapezoidal Rule with n = %d.\n", n);
  fprintf("Estiamted Integral:          %.15f\n", I);
  fprintf("True value:                  %.15f\n", ANS);
  fprintf("True relative percent error: %.15f%%\n\n\n", 100*abs( (ANS-I)/ANS ));

  %// Plots-----------------------------------------
  x_origin = linspace(from-0.5, to+0.5, 100);
  y_origin = f(x_origin);
  x_trap = [from from    x    to    to];
  y_trap = [0    f(from) f(x) f(to) 0 ];

  hold on;
  plot(x_origin, y_origin);
  plot(x_trap, y_trap, '-');
  legend("Original function", "Approximated Trapezoids");
  title("Composite Trapezoidal Rule, n=4");

- Result

  >> HW08_Q1_c
  Composite Trapezoidal Rule with n = 2.
  Estiamted Integral:          16.358608410233252
  True value:                  16.566370614359172
  True relative percent error: 1.254120223205915%
  
  >> HW08_Q1_c
  Composite Trapezoidal Rule with n = 4.
  Estiamted Integral:          16.514833818250274
  True value:                  16.566370614359172
  True relative percent error: 0.311092859797716%

In addition to $n=2$ and $n=4$, I also plot the result with $n=3$, trying to see more from this example.
The method single application is considered as $n=1$ Q1(b) when it comes to composite method. The relative error of Q1(b) is around $5.1816\%$. With four times of $n$ (from $1$ in to $4$), the relative percent error drops to approximately 16 times lower. However, with doubled $n$ (from $1$ in to $4$), it only drops to approximately 16 times lower.
| Full-scope | Zoom-in |
| ------------------------------------ | ------------------------------------ |
|  |  |
|  |  |
|  |  |

### (d) Single application of Simpson’s $1/3$ rule
- Code

  %// Declartions
  f = @(x, varargin) 8 + 4.*cos(x); %// function to integral
  from = 0;                         %// upper limit of integral
  to   = pi/2;                      %// lower limit of integral
  ANS  = (8*to + 4*sin(to)) - (8*from + 4*sin(from));
  m    = (to - from) / 2;

  %// Single application Simpson's 1/3 Rule
  I = m/3 * (f(from) + f(to) + 4*f(m));

  %// Print the result -----------------------------------
  fprintf("Single application Simpson's 1/3 Rule.\n");
  fprintf("Estiamted Integral:          %.15f\n", I);
  fprintf("True value:                  %.15f\n", ANS);
  fprintf("True relative percent error: %.15f%%\n\n\n", 100*abs( (ANS-I)/ANS ));


  %// Plots-----------------------------------------------
  x_pts    = linspace(from-0.5, to+0.5, 100);
  y_origin = f(x_pts);
  a = from;
  b = to;
  m = (a+b) / 2;
  f_lagrange = @(x) f(a).*( (x-m).*(x-b)/((a-m).*(a-b)) ) ...
                  + f(m).*( (x-a).*(x-b)/((m-a).*(m-b)) ) ...
                  + f(b).*( (x-m).*(x-a)/((b-a).*(b-m)) ) ;
  y_trap = f_lagrange(x_pts);

  hold on;
  plot(x_pts, y_origin);
  plot(x_pts, y_trap, '-');
  plot([from from], [0 f(from)], 'g-');
  plot([to to],     [0 f(to)], 'g-');
  legend("Original function", "Lagrange Interpolation", "Integral Limits");
  title("Simpson's 1/3 Rule");

- Result

  >> HW08_Q1_d
  Single application Simpson's 1/3 Rule.
  Estiamted Integral:          16.575490124328013
  True value:                  16.566370614359172
  True relative percent error: 0.055048327609767%

Unlike trapezoidal rule that always under-estiamtes a segement of a concave function, Simpson's 1/3 rule either under or over-esitamtes a segement of a concave function within different intervals.
It evaluates integral of Lagrange interpolation. For this case, it under-estimates at the interval $[0,\frac{\pi}{4})$ and over-estimates at the interval $(\frac{\pi}{4},0]$, decreasing the error of the whole integral range $[0,\frac{\pi}{2}]$.
| Full-scope | Zoom-in |
| ------------------------------------- | ------------------------------------ |
|  |  |

### (e) Composite Simpson’s $1/3$ rule with $n=4$
- Code
> Simpson_13.m is listed in appendix.

  %// Declartions
  f = @(x, varargin) 8 + 4.*cos(x); %// function to integral
  from = 0;                         %// upper limit of integral
  to   = pi/2;                      %// lower limit of integral
  n    = 4;                         %// how many segaments
  ANS  = (8*to + 4*sin(to)) - (8*from + 4*sin(from));

  %// Composite Simpson's 1/3 Rule
  I = Simpson_13(f, from, to, n);

  %// Print the result -----------------------------------
  fprintf("Composite Simpson's 1/3 Rule with n = %d.\n", n);
  fprintf("Estiamted Integral:          %.15f\n", I);
  fprintf("True value:                  %.15f\n", ANS);
  fprintf("True relative percent error: %.15f%%\n\n\n", 100*abs( (ANS-I)/ANS ));

- Result

  >> HW08_Q1_e
  Composite Simpson's 1/3 Rule with n = 4.
  Estiamted Integral:          16.566908954255947
  True value:                  16.566370614359172
  True relative percent error: 0.003249594671677%

### (f) Simpson’s $3/8$ rule
- Code

  %// Declartions
  f = @(x, varargin) 8 + 4.*cos(x); %// function to integral
  from = 0;                         %// upper limit of integral
  to   = pi/2;                      %// lower limit of integral
  ANS  = (8*to + 4*sin(to)) - (8*from + 4*sin(from));
  h    = (to - from) / 3;

  %// Simpson's 3/8 Rule
  I = h*3/8 * (f(from) + f(to) + 3*f((2*from+to)/3) + 3*f((from+2*to)/3))ㄤ

  %// Print the result -----------------------------------
  fprintf("Simpson's 3/8 Rule.\n");
  fprintf("Estiamted Integral:          %.15f\n", I);
  fprintf("True value:                  %.15f\n", ANS);
  fprintf("True relative percent error: %.15f%%\n\n\n", 100*abs( (ANS-I)/ANS ));

- Result

  >> HW08_Q1_f
  Simpson's 3/8 Rule.
  Estiamted Integral:          16.570390307616290
  True value:                  16.566370614359172
  True relative percent error: 0.024264175604238%

### (g) Composite Simpson’s rule with $n=5$.
- Code
> Simpson_38.m is listed in appendix.

  %// Declartions
  f = @(x, varargin) 8 + 4.*cos(x); %// function to integral
  from = 0;                         %// upper limit of integral
  to   = pi/2;                      %// lower limit of integral
  ANS  = (8*to + 4*sin(to)) - (8*from + 4*sin(from));
  n    = 5;
  h    = (to-from) / n;

  %// Mix Simpson's 1/3 and 3/8 Rule to do piecewise integral estimatation

  %// 1/3 Rule
  N = 2;
  a = from;
  b = from + N*h;
  I_31 = Simpson_13(f, a, b, N);

  %// 3/8 Rule
  N = 3;
  a = to - N*h;
  b = to;
  I_38 = Simpson_38(f, a, b, N);

  I = I_31 + I_38;

  %// Print the result -----------------------------------
  fprintf("Mix Simpson's 1/3 and 3/8 Rule with n = %d.\n", n);
  fprintf("Estiamted Integral:          %.15f\n", I);
  fprintf("True value:                  %.15f\n", ANS);
  fprintf("True relative percent error: %.15f%%\n\n\n", 100*abs( (ANS-I)/ANS ));

- Result

  >> HW08_Q1_g
  Mix Simpson's 1/3 and 3/8 Rule with n = 5.
  Estiamted Integral:          16.566704953211904
  True value:                  16.566370614359172
  True relative percent error: 0.002018178033766%

### Error table
Besides the numbers of segmemts $n$ assigned with composite trapezoidal rule and composite Simpson's $1/3$ rule, I also tried some extra $n$.
For composite methods, bigger $n$ yields smaller error in attempts I tried.
> The reason that I say *"in attempts I tried"* is that,
> for big enough of $n$, the error may stop dropping due to the limit of the data type.

The simpson's rule outperforms trapezoidal rule in this case. But I would say that simpson's rule is quite function-shape dependent. For example, estimating integral of a step function, interpolation yields a big error and oscillates with increasing $n$. In the other hand, trapezoidal rule could probabaly yield an a-okay integral estimation with big enough of $n$.
| Method | n | True relative percent error(%) |
|---------------------------|---|--------------------------------|
| Analytically | | 0 |
| Single Trapezoidal | | 5.181625875652982 |
| Composite Trapezoidal | 1 | 5.181625875652982 |
| Composite Trapezoidal | 2 | 1.254120223205915 |
| Composite Trapezoidal | 3 | 0.554168052313233 |
| Composite Trapezoidal | 4 | 0.311092859797716 |
| Single Simpson' s $1/3$ | | 0.055048327609767 |
| Composite Simpson's $1/3$ | 2 | 0.055048327609767 |
| Composite Simpson's $1/3$ | 4 | 0.003249594671677 |
| Composite Simpson's $1/3$ | 6 | 0.000635315271041 |
| Simpson's $3/8$ | | 0.024264175604238 |
| Composite Simpson's | 5 | 0.002018178033766 |

## Q2 Use Romberg integration to evaluate
$$
\int _{\:0}^{2}\frac{e^x\sin \:\left(x\right)}{1+x^2}dx
$$
to an accuracy of $\varepsilon _s=0.5\%$. Your results should be presented in the form of
$O(h^2)\to O(h^4) \to O(h^6) \to O(h^8)$
- Code
The function file `romberg.m` is listed in the appendix.

  clear
  clear all
  close all
  format long

  f = @(x, varargin) exp(x).*sin(x) ./ (1+x.^2);
  a = 0;          %// Integral lower limit
  b = 2;          %// Integral upper limit
  max_iter = 3;   %// romberg of O( h^((max_iter+1)*2) )
  es = 0.09974;   %// Approximation error end condition, unit: percentage

  [I_est,I_,ea,iter] = romberg(f,a,b,es,max_iter,[]);

  %// Print the result
  fprintf("Integrals estimated through romberg:\n");
  fprintf("I_ = \n");
  disp(I_);
  fprintf("The estimated integral of O(h^%d)\n", (iter+1)*2);
  fprintf("I_est = \n");
  disp(I_est);
  fprintf("Es:%g%%, Ea:%g%%\n", es, ea);

- Result
> The approximation error drops $\varepsilon _a$ to $0.997471\%$ at the stage of O($h^6$), while the required order is $O(h^8)$.
> I set $\varepsilon _s=0.09974\%$. It can therefore proceed to $O(h^8)$.

  >> HW08_Q2
  Integrals estimated through romberg:
  I_ = 
     1.343769939485650   1.972826837947778   1.941836048413208   1.939959720716542
     1.815562613332246   1.943772972759119   1.939989038336802                   0
     1.911720382902401   1.940225534238197                   0                   0
     1.933099246404248                   0                   0                   0

  The estimated integral of O(h^8)
  I_est = 
     1.939959720716542

  Es:0.09974%, Ea:0.00151125%

The column of $O(h^2)$ is calculated by trapezoidal rule with $n=1,2,4$ and $8$, where $n$ is the evenly spaced count.
This is for `Romberg` to meet th rule of `Richardon Extrapolation`: $h_{i+1}=h_i/2$
The adavantage of `Romberg` is that it doesn't need too much iterations and yields a better result based on two estimation of different step sizes.
For example, to yield the result of $\text{O} (h^8)$:
　　$I_{est}=1.939959720716542$,
`Romberg` takes
　　$1+2+4+8=15 \text{ iterations to estimate O}(h^2)$
　　$3+2+1=6 \text{ iterations to estimate }O(h^2), O(h^4), O(h^8)$
　　$\text{Total 21 iterations}$,
while trapezoidal rule takes
　　$\text{50 iterations}$
to yield the result that is really close to $I_{est}$

  >> trap(f,a,b,50,[])

  ans =

     1.939950512231325

## Q3 Develop a script to generate the following function in which both independent variables ranging from $-3$ to $3$
### (a) $f(x,y)=e^{-(x^2+y^2)}$
$f(x,y)=0$ as $(x^2+y^2) \to \infty$. This function converges.
- Code

  f = @(x, y) exp( -(x.^2+y.^2) );

  START = -3;  %// start of graph of x and y
  END   = 3;   %// end   of graph of x and y
  DOTS  = 500; %// number of slicing on x, y, z

  x_1D = linspace(START, END, DOTS);
  y_1D = linspace(START, END, DOTS);

  %// [X,Y] = meshgrid([1 2 3],[4 5 6])
  %// X = [1 2 3; 1 2 3; 1 2 3]
  %// Y = [4 4 4; 5 5 5; 6 6 6]
  [x_2D, y_2D] = meshgrid(x_1D, y_1D);
  Z = f(x_2D, y_2D);

  mesh(x_2D, y_2D, Z);
  title('f(x,y) = e^{-(x^2+y^2)}');
  xlabel('x'); ylabel('y'); zlabel('f(x,y)');
  view(45,30); %// change camera view point

- Result
|  |  |
| ------------------------------------ | ------------------------------------ |
|  |  |

### (b) $f(x,y)=xe^{-(x^2+y^2)}$
Since that the power of $e^{-(x^2+y^2)}$ decreases faster than the growing rate of $x$, this function converges as $(x^2+y^2) \to \infty$.
- Code
The code is pretty much same as Q3(a). So only the modified lines are listed, comparing to the code of Q3(a).

  f = @(x, y) x.*exp( -(x.^2+y.^2) );

  title('f(x,y) = xe^{-(x^2+y^2)}');

- Result
|  |  |
| ------------------------------------ | ------------------------------------ |
|  |  |
|  | |

## Q4 Heun’s method
Develope an M-file to solve a single ODE by Heun’s method with iteration. Design the M-file so that it creates a plot of the results.
From chapter 22, we have
　　$y_{i+1}=y_i+h\times \text{Slope}$,
where $\text{Slope}$ in Euler's method is estimated directly from at the beginning of the interval.
　　$\text{Slope}_\text{euler} =\dfrac{dy}{dt}$
In Heun's method, $\text{Slope}$ is determined by the average of the derivatives at beginning and predicted ending of the interval.
　　$\text{Beginning of an interval: }(t_i,y_i)$
　　$\text{Predicted ending of an interval: }(t_i+h,y_i+hf(t_i,y_i))$
　　$\text{Slope}_\text{heun}=\dfrac{1}{2}(f(t_i,y_i)+f(t_i+h,y_i+hf(t_i,y_i)))$,
　　where $f(t,y)=\dfrac{dy}{dt}$
Subsitue $\text{Slope}$ with $\text{Slope}_\text{heun}$ , the formula of Heun's method is acquired:
　　$y_{i+1}=y_i+\dfrac{h}{2}(f(t_i,y_i)+f(t_i+h,y_i+hf(t_i,y_i)))$
- Code
This is modified from `eulode.m`
`heunODE.m`:

  function [t,y] = heunODE(dydt, tspan, y0, h)
  %// heunODE: Heun ODE solver
  %//   [t,y] = eulode(dydt, tspan, y0, h, p1, p2,...)
  %//   `   uses Heun'S method to INTEGRATE an ODE
  %//       (uses the slope at the beginning of the stepsize to graph the
  %//       function.)
  %// Input:
  %//   dydt    = name of hte M-file that evaluates the ODE
  %//   tspan   = [ti,tf] where ti and tf = initial and final values of
  %//               independent variables
  %//   y0      = initial value of dependent variable
  %//   h       = step size
  %//   p1,p2   = additional parameter used by dydt
  %// Output:
  %//   t = vector of independent variable
  %//   y = vector of solution for dependent variable
  if nargin<4, error('at least 4 input arguments required'), end
  ti = tspan(1); tf = tspan(2);
  if ~ (tf>ti), error('upper limit must be greater than lower limit'), end
  t = (ti:h:tf)'; 
  n = length(t);
  %// if necessary, add an additional value of t 
  %// so that range goes from t=ti to tf
  if t(n)<tf
      t(n+1) = tf;
      n = n+1;
      t(n)=tf;
  end
  y = y0*ones(n,1); %// preallocate y to improve efficiency
  for i = 1:n-1 %// implement Heun's Method
      y(i+1) = y(i) + (t(i+1)-t(i))/2 * (...
        dydt(t(i),y(i)) + ...
        dydt(t(i+1), y(i)+h*dydt(t(i),y(i)))...
        );
  end

Call `heunODE.m` to plot. This code is also used in Q5.

  clear
  clear all
  close all
  format long

  %// t_range = [0 2];
  %// y = @(t) exp(t);
  %// dydt = @(t, y) y; %// dy/dt
  %// h = 0.1; %// Steps, bigger step yields bigger error

  %// t_range = [0 6.28];
  %// y = @(t) sin(t);
  %// dydt = @(t, y) cos(t);
  %// h = 0.1; %// Steps, bigger step yields bigger error

  t_range = [-10 10];
  y = @(t) t.*sin(t);
  dydt = @(t, y) sin(t) + t*cos(t);
  h = 0.2; %// Steps, bigger step yields bigger error


  y0 = y(t_range(1));
  %// [t_aprox, y_aprox] = eulode(dydt, t_range, y0, h);
  [t_aprox, y_aprox] = heunODE(dydt, t_range, y0, h);
  %// [t_aprox, y_aprox] = midpointODE(dydt, t_range, y0, h);

  hold on;
  t_plots = linspace(t_range(1), t_range(2), 200);
  plot(t_plots, y(t_plots));
  plot(t_aprox, y_aprox);
  legend("True function", "Approximated", 'location', 'best');
  xlabel("t"); ylabel("y");
  f_str = func2str(dydt);
  tit = sprintf("$y'=%s, t=%d\\sim %d, h=%g$", f_str(7:end), t_range(1), t_range(2), h);
  title(tit,'Interpreter','latex');

- Result
I plot the results and the results of Euler's method with same step sizes are also listed for comparison.
> If there is only one red curve in a graph, it's because the approximated and the true function are really close and the true function is covered.

| Heun's method | Euler's method |
|:------------------------------------:|:------------------------------------:|
|  |  |
|  |  |
|  |  |
Approximated functions of Heun's method are mostly overlapped with true functions in the plots.
It's obvious that Heun's method yields better results than Euler's method with same conditions(initial value, range, step size).

## Q5 Midpoint method
Develop an M-file to solve a single ODE by the midpoint method. Design the M-file so that it creates a plot of the results.
Simlar to $\text{Slope}$ improvement of Euler's method in Heun's, but the derivatives is determine at the middle of an interval rather than at the beginning like Euler's method.
　　$\text{Beginning of an interval: }(t_i,y_i)$
　　$\text{Middle of an interval: }\left( t_i+\dfrac{h}{2},y_i+\dfrac{h}{2}f(t_i,y_i)\right)$
　　$\text{Slope}_\text{mid}=f\left( t_i+\dfrac{h}{2},y_i+\dfrac{h}{2}f(t_i,y_i)\right)$,
　　where $f(t,y)=\dfrac{dy}{dt}$
Subsitue $\text{Slope}$ with $\text{Slope}_\text{mid}$ , the formula of Midpoint method is acquired:
　　$y_{i+1}=y_i+hf\left( t_i+\dfrac{h}{2},y_i+\dfrac{h}{2}f(t_i,y_i)\right)$
- Code
This is modified from `eulode.m`.
`midpointODE.m`

  function [t,y] = midpointODE(dydt, tspan, y0, h)
  %// midpointODE: midpoint ODE solver
  %//   [t,y] = midpointODE(dydt, tspan, y0, h, p1, p2,...)
  %//   `   uses midpoint method to INTEGRATE an ODE
  %//       (uses the slope at the beginning of the stepsize to graph the
  %//       function.)
  %// Input:
  %//   dydt    = name of hte M-file that evaluates the ODE
  %//   tspan   = [ti,tf] where ti and tf = initial and final values of
  %//               independent variables
  %//   y0      = initial value of dependent variable
  %//   h       = step size
  %//   p1,p2   = additional parameter used by dydt
  %// Output:
  %//   t = vector of independent variable
  %//   y = vector of solution for dependent variable
  if nargin<4, error('at least 4 input arguments required'), end
  ti = tspan(1); tf = tspan(2);
  if ~ (tf>ti), error('upper limit must be greater than lower limit'), end
  t = (ti:h:tf)'; 
  n = length(t);
  %// if necessary, add an additional value of t 
  %// so that range goes from t=ti to tf
  if t(n)<tf
      t(n+1) = tf;
      n = n+1;
      t(n)=tf;
  end
  y = y0*ones(n,1); %// preallocate y to improve efficiency
  for i = 1:n-1 %// implement midpoint Method
      y(i+1) = y(i) + (t(i+1)-t(i)) * dydt(...
        t(i) + (t(i+1)-t(i))/2, ...
        y(i) + (t(i+1)-t(i))/2 * dydt(t(i), y(i))...
        );
  end

The part of calling `midpoint.m` is almost the same as the one in Q5. The only difference is that `heunODE()` is not called but `midpoint.m`.

  %// [t_aprox, y_aprox] = eulode(dydt, t_range, y0, h);
  %// [t_aprox, y_aprox] = heunODE(dydt, t_range, y0, h);
  [t_aprox, y_aprox] = midpointODE(dydt, t_range, y0, h);

- Result

| Midpoint method | Euler's method |
|:------------------------------------:|:------------------------------------:|
|  |  |
|  |  |
|  |  |
It's obvious that midpoint method yields a result that is more accurate than Euler's method.

## Q6 Given dy/dt:
$\dfrac{dy}{dt} = -100000y+99999e^{-t}$
It can be solved analytically.
Use the form
　　$y'(t)+p(t)y(t)=q(t)$
　　$y'+10^5y=99999e^{-t}$,　$p(t)=10^5$,　$q(t)=99999e^{-t}$
Find a $u(t)$ that satiesfies
　　$(u(t)y(t))'=u(t)q(t)$
Solve $u(t)$
$$
\begin{aligned}
(uy)' &= u'y+uy' \\
&= u'y+u(-10^5y+99999e^{-t}) \\
&= uq \\
&= u\times 999999e^{-t} \\
\end{aligned}
$$
$$\begin{aligned}
&u'y+u(-10^5y+99999e^{-t})=u\times 999999e^{-t} \\
&u'y+u(-10^5y)=0 \\
&u'y=10^5uy \\
&u'=10^5u, \text{ For } y\neq 0 \\
&u=u(t)=e^{10^5t}, \text{ For } u(0)=1
\end{aligned}
$$
　　$(e^{10^5t}y)'=e^{10^5t}\times99999e^{-t}=99999e^{99999t}$
　　$y=e^{-t}+ce^{-10^5t}$, where $c$ is a constant determined by initial value.
### (a) Find step size
Estimate the step size required to maintain stability using the explicit Euler method.
The explicit Euler method:
$$y_{i+1}=y_i+\dfrac{dy}{dt}(t_i, y_i)\times h,$$
where $h$ is the step size.
To evaluate convergence,
$$\begin{aligned}
&\text{let } t\to \infty \\
&\dfrac{dy}{dt}=-100000y+99999e^{-t}=-10^5y
\end{aligned}$$
Subsitute $\dfrac{dy}{dt}\bigg|_{t\to\infty}$ with explicit Euler method
$$
\begin{aligned}
y_{i+1}&=y_i-10^5y_i\times h\\
y_{i+1}&=y_i(1-10^5h)
\end{aligned}
$$
For $y_i$ to converge when $t\to \infty$, the change of $y_i$ should be no bigger then $1$.
That is,
$$\begin{aligned}
|1 &- 10^5h|<1 \\
-1<1 &- 10^5h<1 \\
-2<0 &- 10^5h<0 \\
&h<\dfrac{2}{10^5}
\end{aligned}$$
To maintain the stability using the explicit Euler method, the step size $h$ should be smaller than $\dfrac{2}{10^5}$
- Code

clear
clear all
close all
format long

t_range = [0 1];
c = -1; %// constant
y = @(t) exp(-t) + c .* exp(-1e5.*t);
dydt = @(t, y) -1e5.*y + 99999.*exp(-t);
h = 1*1e-5; %// Steps, bigger step yields bigger error

y0 = y(t_range(1));
[t_aprox, y_aprox] = eulode(dydt, t_range, y0, h);

hold on;
plot(t_aprox, y(t_aprox), 'g');
plot(t_aprox, y_aprox);
legend("True", "Approximated", 'location', 'best');

xlabel("t"); ylabel("y");
tit = sprintf("$h=%g$", h);
title(tit,'Interpreter','latex');

- Result
| Full-scope | Zoom-in |
|--------------------------------------|--------------------------------------|
|  |  |
|  |  |
|  |  |
I plot both approximated function and original function solve with online ODE solver.
> I tried to solve it by myself, but I failed with hours of efforts. And when I finally figured it out, I don't have time to include the solution in this homework.

The result shows that for step size $h\geq2\times10^{-5}$, explicit euler method failed to converge and ends up with either oscillation or $\pm\infty$.
In the other hand, for step size $h\leq2\times10-5$, such as $h=1.99\times 10^{-5}$, the result oscillates in the beginning but matches the true function pretty well.

### (b) Initial value
If $y(0)=0$, use the implicit Euler to obtain a solution from $t=0$ to $t=2$ using step size of $0.1$.
Subsitute
　　$y(0)=0$　to　$y=e^{-t}+ce^{-10^5t}$,
　　$c=-1$
　　$y=e^{-t}-e^{-10^5t}$
Formula of implicit euler method:
　　$y_{i+1}=y_i+f(t_{i+1}, y_{i+1})h$,
where $f(t_{i+1}, y_{i+1})=\dfrac{dy}{dt}$ and it is known.
Subsitute it with $-100000y+99999e^{-t}$
$$\begin{aligned}
y_{i+1} &=y_i+(-10^5y_{i+1}+99999e^{-t_{i+1}})h \\
(1+10^5h)y_{i+1}&=y_i+e^{-t_{i+1}}h \\
y_{i+1} &= \dfrac{y_i+e^{-t_{i+1}}}{1+10^5h}
\end{aligned}
$$
The formula of implicit euler method for the given $dy/dt$ is acquired.
The corresponding code in iterration:

y_est(i+1) = ( y_est(i)+99999*exp(-t_span(i+1))*h ) / (1+1e5*h);

- Code

clear
clear all
close all
format long

t_range = [0 2];
c = -1; %// constant
y = @(t) exp(-t) + c .* exp(-1e5.*t);
dydt = @(t, y) -1e5.*y + 99999.*exp(-t);
h = 0.1; %// Steps, bigger step yields bigger error

%// Use implicit Euler ODE
t_span = t_range(1):h:t_range(2);
n = length(t_span);
y_est = zeros(1, n);
y_est(1) = y(t_range(1)); %// t_range(1) is beginging of t_span
for i = 1:(n-1)
  y_est(i+1) = ( y_est(i)+99999*exp(-t_span(i+1))*h ) / (1+1e5*h);
end

hold on;
plot(t_span, y_est, 'b');
t_plots = linspace(t_range(1), t_range(2), 200);
plot(t_plots, y(t_plots), 'g');
legend("Approximated", "True", 'location', 'best');

- Result
| | |
| ------------------------------------ | ------------------------------------ |
|  |  |
Unlike the explicit one, this implicit method converges with a high tolerance step size.
But it's harder to solve an ODE with implicit euler mothod than the explicit one automatically.

## Appendix
### `trap.m`

function I = trap(func,a,b,n,varargin)
%// trap: composite trapezoidal rule quadrature
%//   I=trap(func,a,b,n,p1,p2,...):
%//                composite trapezoidal rule
%// input:
%//   func=name of fuction to be integrated
%//   a, b=integration limits
%//   n=number of segments (default=100)
%//   p1, p2,...=additional parameters used by function
%// output:
%//   I=integral estimate
if nargin<3,error('at least 3 input arguments required'),end
if ~(b>a),error('upper bound must be greater than lower'),end
if nargin<4||isempty(n),n=100;end
x=a; h=(b-a)/n;
s=func(a,varargin{:});
for i=1:n-1
    x=x+h;
    s=s+2*func(x,varargin{:});
end
s=s+func(b,varargin{:});
I=(b-a)*s/(2*n);

### `romberg.m`

function [q,I_,ea,iter]=romberg(func,a,b,es,maxit,varargin)
%// romberg: Romberg integration quadrature
%// q = romberg(func,a,b,es,maxit,p1,p2,...):
%// Romberg integration.
%// input:
%//   func = name of function to be integrated
%//   a, b = integration limits
%//   es = desired relative error (default = 0.000001%)
%//   maxit = maximum allowable iterations (default = 30)
%//   pl,p2,... = additional parameters used by func
%// output:
%//   q = integral estimate
%//   I_ = integrals through the iteration
%//   ea = approximate relative error (%)
%//   iter = number of iterations
  if nargin<3,error('at least 3 input arguments required'),end
  if nargin<4|isempty(es), es=0.000001;end
  if nargin<5|isempty(maxit), maxit=30;end
  n = 1;
  I(1,1) = trap(func,a,b,n,varargin{:});
  iter = 0;
  while iter<maxit
    iter = iter+1;
    n = 2^iter;
    I(iter+1,1) = trap(func,a,b,n,varargin{:});
    for k = 2:iter+1
      j = 2+iter-k;
      I(j,k) = (4^(k-1)*I(j+1,k-1)-I(j,k-1))/(4^(k-1)-1);
    end
    ea = abs((I(1,iter+1)-I(2,iter))/I(1,iter+1))*100;
    if ea<=es, break; end
  end
  I_ = I;
  q = I(1,iter+1);
end

### `Simpson_13.m`

function I = Simpson_13(f, a, b, n)
%// Usage:
%//   I = Simpson_31(func,a,b,n,p1,p2,...):
%// input:
%//   f    = name of fuction to be integrated
%//   a, b = integration limits
%//   n    = number of segments, should be multiple of 3
%// output:
%//   I=integral estimate

if rem(n,2)==1
   error('\n Argument n invalid. n should be multiple even.'); 
end

%// Composite Simpson's 1/3 Rule
h      = (b-a) / n; %// aka h
x      = linspace(a+h, b-h, n-1);
x_odd  = x(1:2:end);    %// x(i) for all mod(i,2)==1
x_even = x(2:2:end);    %// x(i) for all mod(i,2)==0 && 2<=i<=n-2

I = h/3 * (f(a) + f(b) + 4*sum(f(x_odd)) + 2*sum(f(x_even)));

### `Simpson_38.m`

function I = Simpson_38(f, a, b, n)
%// Ref:
%//   https://www.mathworks.com/matlabcentral/answers/587858-3-8-simpson-s-rule
%// Usage:
%//   I = Simpson_38(func,a,b,n):
%// input:
%//   f    = name of fuction to be integrated
%//   a, b = integration limits
%//   n    = number of segments, should be multiple of 3
%// output:
%//   I=integral estimate

if rem(n,3)~=0
   error('\n Argument n invalid. n should be multiple of 3'); 
end
x = linspace(a,b,n+1);
I = 3*(b-a)/8/n*sum(f(x).*[1,3,3,repmat([2,3,3],1,(n-3)/3),1]);

### `eulode.m`

function [t,y] = eulode(dydt, tspan, y0, h)
%eulode: Euler ODE solver
%//   [t,y] = eulode(dydt, tspan, y0, h, p1, p2,...)
%//   `   uses EULER'S method to INTEGRATE an ODE
%//       (uses the slope at the beginning of the stepsize to graph the
%//       function.)
%Input:
%//   dydt    = name of hte M-file that evaluates the ODE
%//   tspan   = [ti,tf] where ti and tf = initial and final values of
%//               independent variables
%//   y0      = initial value of dependent variable
%//   h       = step size
%//   p1,p2   = additional parameter used by dydt
%Output:
%//   t = vector of independent variable
%//   y = vector of solution for dependent variable
if nargin<4, error('at least 4 input arguments required'), end
ti = tspan(1); tf = tspan(2);
if ~ (tf>ti), error('upper limit must be greater than lower limit'), end
t = (ti:h:tf)'; 
n = length(t);
%// if necessary, add an additional value of t 
%// so that range goes from t=ti to tf
if t(n)<tf
    t(n+1) = tf;
    n = n+1;
    t(n)=tf;
end
y = y0*ones(n,1); %// preallocate y to improve efficiency
for i = 1:n-1 %// implement Euler's Method
    y(i+1) = y(i) + dydt(t(i),y(i))*(t(i+1)-t(i));
end