---
# System prepended metadata

title: Numerical Methods Homework-6
tags: [' Numerical Methods', homeworks]

---

---
tags: homeworks, Numerical Methods
---
# Numerical Methods Homework-6
## B10702026 林琨霖
[Online link of this homework](https://hackmd.io/iZg_1Jw8T4W0p0f_MlSxFw)
It's recommanded to review this homework through the link above. Then content is same but with a better layout, clickable index and syntax highlighting.
> For the online link, you can click `...` at the upper right corner to check for the last edition time, making sure this homework is finished before the dead line.
> ![](https://i.imgur.com/Liwq2th.png =600x)

## Q1 Directional derivative
Find the directional derivative of:
$$
f(x,y)=2x^2+y^2
$$
at $x=2$ and $y=4$ in the direction of $h=3i+2j$.
### Calculation
- Declartions
  ```c
  % //symbolic declartion
  syms x;
  syms y;
  f = 2.*x.^2 + y.^2;

  % //Direction: h=3i+2j
  h = [3 2];

  % //At (p_x, p_y)
  p_x = 2;
  p_y = 4;
  ```
- Gradient:
  $\frac{\partial f}{\partial x}=4x$
  $\frac{\partial f}{\partial y}=2y$
  $\nabla f=(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y})$
  $\nabla f(2,4)=(8,8)=8i+8j$
  ```c
  partial_x = diff(f, x);
  partial_y = diff(f, y);
  grad = [partial_x partial_y];
  ```
- Directional derivative
  $D(f)_{h}=\nabla f \cdot \frac{h}{\big| h\big| }$
  $D(f)_{h}=(8,8)\cdot \frac{(3,2)}{\big| (3,2)\big| }$
  $D(f)_{h}=(8,8)\cdot \frac{(3,2)}{\sqrt {13}}$
  $D(f)_{h}=\frac{40}{\sqrt {13}}=11.09400392$
  ```c
  x = p_x;
  y = p_y;
  answer = dot(subs(grad), h)/ sqrt(sum(h.^2))
  double(answer) % //Evaluate symbolic expression
  ```
### Code
  ```c=
  clear
  clear all
  close all
  format long


  % //symbolic declartion
  syms x;
  syms y;
  f = 2.*x.^2 + y.^2;

  % //Direction: h=3i+2j
  h = [3 2];

  % //At (p_x, p_y)
  p_x = 2;
  p_y = 4;

  % //Graph tables
  graph_x = linspace(p_x-2, p_x+2, 100);
  graph_y = linspace(p_y-2, p_y+2, 100);

  % //Gradient of f
  partial_x = diff(f, x);
  partial_y = diff(f, y);
  grad = [partial_x partial_y];

  % //Evaluate directional derivative
  x = p_x;
  y = p_y;
  answer = dot(subs(grad), h)/ sqrt(sum(h.^2))
  double(answer) % //Evaluate symbolic expression


  % //3D plots -----------------------------------
  hold on;
  title('D(f)_{h}, (x,y)=(2,4), h=3i+2j');
  xlabel('x-axis');
  ylabel('y-axis');
  zlabel('f(x,y)');
  view(45,45);
  pt = plot3(x, y, double(subs(f)), 'or');


  [x_2D, y_2D] = meshgrid(graph_x, graph_y);
  x = x_2D;
  y = y_2D;
  mesh(x_2D, y_2D, double(subs(f)));

  x = p_x; y = p_y;
  % //quiver3(X, Y, Z, dirX, dirY, dirZ, arrow_length);
  quiver3(x, y, double(subs(f)), ...
          h(1)/sqrt(sum(h.^2)), h(2)/ sqrt(sum(h.^2)), 0, 1);

  legend(pt(1), 'f(2, 4)', 'Location', 'best');
  legend('mesh of f(x,y)');
  legend('Unit vector of h');
  ```
### Result
  ```
  >> HW06_Q1
  answer =
    (40*13^(1/2))/13

  ans =
    11.094003924504582
  ```
  The directional derivative of $f(x,y)$ at $x=2$ and $y=4$ in the direction of $h=3i+2j$ is ==11.094003924504582==.
  |                                                                               |                                                                             |
  |-------------------------------------------------------------------------------|-----------------------------------------------------------------------------|
  | x-y-f()![](https://i.imgur.com/QlpHrbj.png =540x300)                             | x-y-f()![](https://i.imgur.com/AzKw7yR.png =540x300)                           |
  | x-f()<br>![](https://i.imgur.com/ctwNkMG.png =540x300)                               | y-f()<br>![](https://i.imgur.com/1zyU576.png =540x300)                             |
  | x-y (Origin at bottom-right corner)<br>![](https://i.imgur.com/i9bcUUI.png =540x300) | x-y (Origin at upper-left corner)<br>![](https://i.imgur.com/N5L3c8Y.png =540x300) |



---
## Q2 Maximize $f(x,y)$
Given:
$$f(x,y)=2.25xy+1.75y-1.5x^2-2y^2$$
Construct and solve a system of linear algebraic equations that maximizes $f(x,y)$. Notice that this is done by setting the partial derivative of $f$ with respect to both $x$ and $y$ to zero.
### Calculation
- Declartions
  ```c
  syms x;
  syms y;
  f = 2.25.*x.*y + 1.75.*y - 1.5*x.^2 - 2.*y.^2;
  ```
- Partial derivative
  $\begin{cases}
    \dfrac{\partial f}{\partial x}=2.25y-3x\\
    \\
    \dfrac{\partial f}{\partial y}=2.25x+1.75-4y
  \end{cases}$
  ```c
  partial_x = diff(f, x);
  partial_y = diff(f, y);
  ```
  Let partials to be zero, we have these two equation
  $\begin{cases}
    -3x+2.25y=0\\
    \\
    2.25x-4y=-1.75
  \end{cases}$
Solve for $(x,y)$, where the extreme value at.
- Find where extreme occurs
  Solving the equation using matrix $Ax=b$. I will replace $x$ with $k$ in case of variable names confusing.
  $$
  A=\begin{bmatrix}-3&2.25\\ 2.25&-4\end{bmatrix},\:k=\begin{bmatrix}x\\ y\end{bmatrix},b=\begin{bmatrix}0\\ -1.75\end{bmatrix}\\
  $$
  ```c
  % //First row of A
  x=1; y=0; A(1,1) = subs(partial_x) - partial_x_const;
  x=0; y=1; A(1,2) = subs(partial_x) - partial_x_const;
  % //Second row of A
  x=1; y=0; A(2,1) = subs(partial_y) - partial_y_const;
  x=0; y=1; A(2,2) = subs(partial_y) - partial_y_const;
  b = [-partial_x_const; -partial_y_const]; % //negative because it's at the opposite of the =
  ```
  Solve for $k$
  $$
  Ak=b, \ k=A^{-1}*b
  $$
  ```c
  k = A\b
  ```
  We have the point where extreme point occurs
  $$
  k=\begin{bmatrix}x\\ y\end{bmatrix}=\begin{bmatrix}\dfrac{21}{37}\\ \dfrac{28}{37}
  \end{bmatrix}
  $$
- Evaluate the extreme value
  $(x,y)=\left( \dfrac{21}{37}, \dfrac{28}{37} \right)$
  $\begin{aligned}
    f\left(\dfrac{21}{37}, \dfrac{28}{37}\right) &=2.25*\dfrac{21}{37}*\dfrac{28}{37}+1.75*\dfrac{28}{37}-1.5*\left(\dfrac{21}{37}\right)^2-2*\left( \dfrac{28}{37}\right) ^2 \\
    &=\dfrac{49}{74}
  \end{aligned}$
  ```c
  x = k(1);
  y = k(2);
  extreme = subs(f)
  ```
### Code
  ```c=
  clear
  clear all
  close all
  format long

  % //symbolic declartion
  syms x;
  syms y;
  f = 2.25.*x.*y + 1.75.*y - 1.5*x.^2 - 2.*y.^2;
  partial_x = diff(f, x);
  partial_y = diff(f, y);

  x = 0; y = 0;
  partial_x_const = subs(partial_x); % //constant of partial_x
  partial_y_const = subs(partial_y); % //constant of partial_y

  % //Ak=b, solve for k, where k=[x;y]
  % //First row of A
  x=1; y=0; A(1,1) = subs(partial_x) - partial_x_const;
  x=0; y=1; A(1,2) = subs(partial_x) - partial_x_const;
  % //Second row of A
  x=1; y=0; A(2,1) = subs(partial_y) - partial_y_const;
  x=0; y=1; A(2,2) = subs(partial_y) - partial_y_const;
  b = [-partial_x_const; -partial_y_const]; % //negative because it's at the opposite of the =
  k = A\b; % //x<--k(1), y<--k(2)

  % //Evluate the extreme value
  x = k(1)
  y = k(2)
  extreme = subs(f)
  extreme = double(extreme)


  % //3D plots -----------------------------------
  hold on;
  title('f(x,y) has extreme valur at (21/37, 28/37)');
  xlabel('x-axis');
  ylabel('y-axis');
  zlabel('f(x,y)');
  view(45,45);
  pt = plot3(x, y, double(subs(f)), 'or');

  % //Graph tables
  graph_x = linspace(subs(x)-2, subs(x)+2, 100);
  graph_y = linspace(subs(y)-2, subs(y)+2, 100);

  [x_2D, y_2D] = meshgrid(graph_x, graph_y);
  x = x_2D;
  y = y_2D;
  pt_mesh = mesh(double(x_2D), double(y_2D), double(subs(f)));

  % //Graph 2D lines
  L1 = @(x) (b(1)-A(1,1).*x) / A(1,2);
  L2 = @(x) (b(2)-A(2,1).*x) / A(2,2);
  x_s = linspace(-1.5, 2, 100);
  pt_L1 = plot(x_s, L1(x_s), '-g');
  pt_L2 = plot(x_s, L2(x_s), '-b');

  legend([pt(1) pt_L1 pt_L2], 'Extreme point', 'Partial x=0', 'Partial y=0', 'Location', 'best');

  ```
### Result
  ```
  >> HW06_Q2
 
  x =
    21/37

  y =
    28/37

  extreme =
    49/74

  extreme =
     0.662162162162162
  ```
  The maximal value of $f(x,y)$ is ==$\dfrac{49}{74}=0.6\overline{62162}$== at ==$x=\dfrac{21}{37}=0.\overline{567}$== and ==$y=\dfrac{28}{37}=0.\overline{756}$==.
  Note that the green and blue lines are at drawn at 2D plane, $z=0$, which is irrelevant. And the intersection $(x,y)$ of these two lines is where the exterme value of $f(x,y)$ occurs.  
  | y-x<br>![](https://i.imgur.com/6CaUmMP.png =430x300)      | x-y<br>![](https://i.imgur.com/Etbre0E.png =430x300)      |
  | --------------------------------------------------------- | --------------------------------------------------------- |
  | x-f(x,y)<br>![](https://i.imgur.com/Uh8qXFu.png =430x300) | y-f(x,y)<br>![](https://i.imgur.com/8Ihlb5M.png =430x300) |
  | ![](https://i.imgur.com/as1mDNa.png =430x300)             | ![](https://i.imgur.com/bDc5eQF.png =430x300)             |
  | ![](https://i.imgur.com/aKtE1dh.png =430x300)             | ![](https://i.imgur.com/3uMJr5Z.png =430x300)             |

---
## Q3 Eigenvalues and Hooke's law
Consider a system with 2 frictionless masses connected to walls by 3 linear elastic springs, where $M_1=5\text{kg}$, $M_2=10\text{kg}$, $k_a=0.1\text{N/m}$, $k_b=11\text{N/m}$, and $k_c=0.2\text{N/m}$.
Plot the positions of both carts by modeling it as a 2*2 matrix and then find the roots of the characteristic polynomial.
![](https://i.imgur.com/8zIYUKg.png)
### Calculation
- Declartion
  The following two equations(corresponding to two masses) are based on *Newton's 2nd law* and *Hooke's law*.
  $\ 　　F=ma=-kx$,
  $\ 　　ma+kx=0$, 
  where $k$ is $\text{spring constant}$ and $x$ is $\text{displacement}$ along x-axis.
  For the two masses, we have  two equation:
  $\ 　　\text{eq1}: m_1a_1+x_1(k_a+k_b)-x_2k_b =0$
  $\ 　　\text{eq2}: m_2a_2-x_1k_b+x_2(k_b+k_c) =0$
  We also know that the a mass with a spring results in sinusoid motion.
  $\ 　　x_i=X_i*\text{sin}(\omega t)$
  $\ 　　a_i=\dfrac{d^2x_i}{dt^2}=-X_i \omega^2 \text{sin}(\omega t)$, where $i=1,2$.
  Subsitute back to $\text{eq1}$ and $\text{eq2}$:
  $\begin{cases}
    -X_1m_1w^2\text{sin}(tw) + X_1\text{sin}(tw)(k_a+k_b) - X_2k_b\text{sin}(tw)=0\\
    \\
    -X_1k_b\text{sin}(tw) + X_2\text{sin}(tw)(k_b+k_c) - X_2m_2w^2\text{sin}(tw)=0
  \end{cases}$
  <br>
  $\begin{cases}
    X_1\text{sin}(tw)(k_a+k_b - m_1w^2) - X_2k_b\text{sin}(tw)=0\\
    \\
    -X_1k_b\text{sin}(tw) + X_2\text{sin}(tw)(k_b+k_c-m_2w^2)=0
\end{cases}$
  ```c
  syms X1; syms X2;
  syms w; syms t;
  syms m1; syms m2;
  syms ka; syms kb; syms kc;

  %// Sin
  x1 = X1.*sin(w.*t); a1 = diff(x1, t, 2);
  x2 = X2.*sin(w.*t); a2 = diff(x2, t, 2);

  %// F = m*a = -k*x
  %// m*a + k*x = 0
  eq1 = m1*a1 + x1*(ka+kb) - x2*(kb)    % equaltion 1
  eq2 = m2*a2 - x1*(kb)    + x2*(kb+kc) % equaltion 2
  ```
  Collect the coefficients of $X_1$ and $X_2$ and form a $2\text{x}2$ matrix:
  $\ 　　Al=\begin{bmatrix}\text{sin}\left(tw\right)\left(k_a+k_b\:-\:m_1w^2\right)&-k_b\text{sin}\left(tw\right)\\ -k_b\text{sin}\left(tw\right)&\text{sin}\left(tw\right)\left(k_b+k_c-m_2w^2\right)\end{bmatrix}$
  <br>
  $\ 　　Al'=\begin{bmatrix}\:\frac{k_a+k_b}{\:m_1}-\:w^2&\frac{-k_b}{m_1}\\ \:\frac{-k_b}{m_2}&\frac{k_b+k_c}{m_2}-w^2\end{bmatrix}$
  <br>
  $\ 　　Al'=A-\lambda I,\ \ I=\begin{bmatrix}1&0\\ 0&1\end{bmatrix},\ \ \lambda=w^2$
 ```c
 %// (A-lambda*I) * [X1; X2] = 0
 %// extract coefficients from eq1 and eq2 to matrix Al
 X1 = 1; X2 = 0;
 Al(1, 1) = subs(eq1);
 Al(2, 1) = subs(eq2);

 X1 = 0; X2 = 1;
 Al(1, 2) = subs(eq1);
 Al(2, 2) = subs(eq2);
 Al
 ```
- Solve
  $\text{Let } \text{det}(Al)=0$
  We have
  $\ 　　\text{sin}(tw)\left( k_ak_b + k_ak_c + k_bk_c + w^2(-k_am_2-k_bm_1-k_bm_2-k_cm_1) + w^4m_1m_2 \right)=0$
  $\ 　　k_ak_b + k_ak_c + k_bk_c + w^2(-k_am_2-k_bm_1-k_bm_2-k_cm_1) + w^4m_1m_2=0$
  Subsitute the given parameters:
  $\ 　　8/25 - 17w^2 + 50w^4 = 0$
  Solve for $w$ then we have the possible angular frequencies for the system:
  $\ 　　w=\dfrac{2\sqrt 2}{5},\ \ \dfrac{\sqrt 2}{10}$, negative frequency not applicable.
  ```c
  %// Solve for omega(which is w)
  eq_w = det(Al);         %// equation of omega
  w_s = solve(eq_w, w);   %// w_s stands for multi roots of omega
  w_s = w_s( w_s>=0 );    %// Exclude the negatives
  ```
### Code
  ```c=
  clear
  clear all
  close all
  format long

  m1=5; m2=10;
  ka=0.1; kb=1; kc=0.2;

  syms X1; syms X2;
  syms w; syms t;

  %// Sin
  x1 = X1.*sin(w.*t); a1 = diff(x1, t, 2);
  x2 = X2.*sin(w.*t); a2 = diff(x2, t, 2);

  %// F = m*a = -k*x
  %// m*a + k*x = 0
  eq1 = m1*a1 + x1*(ka+kb) - x2*(kb)    %// equaltion 1
  eq2 = m2*a2 - x1*(kb)    + x2*(kb+kc) %// equaltion 2

  %// (A-lambda*I) * [X1; X2] = 0
  %// Al = A-lambda*I
  %// extract coefficients from eq1 and eq2 to matrix Al
  X1 = 1; X2 = 0;
  Al(1, 1) = subs(eq1);
  Al(2, 1) = subs(eq2);

  X1 = 0; X2 = 1;
  Al(1, 2) = subs(eq1);
  Al(2, 2) = subs(eq2);
  Al

  %// Solve for omega(which is w and w^2 is lambda)
  eq_w = det(Al);         %// equation of omega
  w_s = solve(eq_w, w);   %// w_s stands for multi roots of omega
  w_s = w_s( w_s>0 );     %// Exclude non-positives
  eig_value = [w_s(1)^2 0; 0 w_s(2)^2];

  %// Solve for eigen factor with omega 1
  w = w_s(1);
  t = 1; %// t=1 for evaluating Al
  A_e = subs( Al );
  v(1,1) = -A_e(2,2) / sqrt( A_e(2,1).^2+A_e(2,2).^2 );
  v(2,1) =  A_e(2,1) / sqrt( A_e(2,1).^2+A_e(2,2).^2 );
  v = double(v);

  %// Solve for eigen factor with omega 2
  w = w_s(2);
  t = 1; %// t=1 for evaluating Al
  A_e = subs( Al );
  v(1,2) = -A_e(1,2) / sqrt( A_e(1,1).^2+A_e(1,2).^2 );
  v(2,2) =  A_e(1,1) / sqrt( A_e(1,1).^2+A_e(1,2).^2 );
  v = double(subs(v));

  %// Show the result
  fprintf("omega^2 =\n");
  fprintf("%17s = %.10f\n", w_s(1)^2, double( w_s(1)^2 ));
  fprintf("%17s = %.10f\n", w_s(2)^2, double( w_s(2)^2 ));
  fprintf("\n\n");
  fprintf("omega =\n");
  fprintf("%17s = %.10f\n", w_s(1), double( w_s(1) ));
  fprintf("%17s = %.10f\n", w_s(2), double( w_s(2) ));
  v

  %// Plots ----------------------------------------------
  for i = [1 2]
    t = linspace(0, 30*i, 150);
    figure(i);
    hold on;
    w = double(w_s(i));
    title("omega^2 = " + num2str(w^2));
    X1 = v(1,1)/v(1,1);
    X2 = v(2,1)/v(1,1);
    plot(t, subs(x1));
    plot(t, subs(x2));
    legend('mass1', 'mass2', 'location', 'best');
    xlabel('t (time)'); ylabel('x (displacement)');
  end
  ```
### Result
  ```
  >> HW06_Q3
 
  eq1 =
      (11*X1*sin(t*w))/10 - X2*sin(t*w) - 5*X1*w^2*sin(t*w)


  eq2 =
      (6*X2*sin(t*w))/5 - X1*sin(t*w) - 10*X2*w^2*sin(t*w)


  Al =
    [(11*sin(t*w))/10 - 5*w^2*sin(t*w),                        -sin(t*w)]
    [                        -sin(t*w), (6*sin(t*w))/5 - 10*w^2*sin(t*w)]

  omega^2 =
               8/25 = 0.3200000000
               1/50 = 0.0200000000

  omega =
      (2*2^(1/2))/5 = 0.5656854249
         2^(1/2)/10 = 0.1414213562

  v =
     0.894427190999916   0.707106781186548
    -0.447213595499958   0.707106781186548
  ```
  The roots of the characteristic polynomial is `omega^2`= ==0.32== or ==0.02==.
  |                              |              $w^2=0.32$              |                         $w^2=0.02$                        |
  |:----------------------------:|:------------------------------------:|:---------------------------------------------------------:|
  | Displacement<br>v.s.<br>Time | ![](https://i.imgur.com/0cTfxR1.png) | Two lines overlapped ![](https://i.imgur.com/4YTbTwW.png) |
  |   Location<br>v.s.<br>Time   | ![](https://i.imgur.com/c5TFLwU.png) |            ![](https://i.imgur.com/t6O8Xbg.png)           |

---
## Q4 Non-linear Equation with Linear Regression 
An inverse investigator has reported the data tabulated below. It is known that such
data can be modeled by the following equation:
$\ 　　x=e^{\left(y-b\right)/a}$
where a and b are parameters. Use a transformation to linearize this equation and
then employ linear regression to determine a and b. Based on your analysis predict $y$ at $x=2.6$.
$\ 　　\begin{array} {|r|r|}\hline x & 1 & 2 & 3 & 4 & 5 \\ \hline y & 0.5 & 2 & 2.9 & 3.5 & 4 \\ \hline  \end{array}$
### Calculation
- Declartion
  ```c
  true_x = [1 2 3 4 5];
  true_y = [0.5 2 2.9 3.5 4];
  n = length(true_x);
  ```
- Transformation
  We have the equation of $x$:
  $\ 　　x=e^{\left(y-b\right)/a}$
  $\ 　　ln(x)=\left(y-b\right)/a$
  $\ 　　y=a*ln(x)+b,\text{ let }M=ln(x)$
  $\ 　　y=a*M+b$
- Regression by formula
  Now it's considered as a linear equation and the formula of linear regresion can be applied.
  $\ 　　a=\dfrac{n\sum{M_i y_i}-\sum{M_i}*\sum{y_i}}{n\sum{M_i^2}-(\sum{M_i})^2}$
  $\ 　　b=y_{avg}-aM_{avg}$
  ```c
  a = ( n*sum(log(true_x).*true_y)-sum(log(true_x))*sum(true_y) ) / ...
       ( n*sum(log(true_x).^2)-(sum(log(true_x))^2 ) );
  b = mean(true_y) - a*mean(log(true_x));
  ```
- Pridiction of $y$ at $x=2.6$
  Once $a$ and $b$ are acquired, we can $y$ evaluate approximately:
  $\ 　　y=a*M+b=a*ln(x)+b$
  $\ 　　y(2.6)=a*ln(2.6)+b$
  ```c
  y_approxi = f_transed(2.6)
  ```
### Code
  ```c=
  clear all
  close all
  format long

  eval_x = 2.6;
  true_x = [1 2 3 4 5];
  true_y = [0.5 2 2.9 3.5 4];
  n = length(true_x);

  % Estiated line: y = a1*M + a0, where M is ln(x)
  a = ( n*sum(log(true_x).*true_y)-sum(log(true_x))*sum(true_y) ) / ...
       ( n*sum(log(true_x).^2)-(sum(log(true_x))^2 ) );
  b = mean(true_y) - a*mean(log(true_x));

  % Before transform: x = e^((y-b)/a)
  % ln(x) = (y-b)/a, y = a*ln(x)+b
  f_transed = @(x) a*log(x) + b;

  % Original funxtion
  fprintf("a=%.15f\n", a);
  fprintf("b=%.15f\n", b);
  y = @(x) a*log(x)+b
  x = @(y) log((y-b)/a)

  approxi_y = f_transed(eval_x)

  % //Plot Of 'Transformed to ln(x)-y'
  figure(1); hold on;
  title('Transformed: y-ln(x)');
  xlabel('ln(x) axis');
  ylabel('y axis');
  plot(log(true_x), true_y, 'o');
  ln_x = linspace(min(log(true_x))-1, max(log(true_x))+1, 100);
  plot(ln_x, a*ln_x + b);
  plot(log(eval_x), approxi_y, '^');
  legend('Given Data', 'Estimated curve', 'x=2.6', 'location', 'best');

  % //Plot Of 'Original form x-y'
  figure(2); hold on;
  title('Original form: x-y');
  xlabel('x axis');
  ylabel('y axis');
  plot(true_x, true_y, 'o');
  x = linspace(min(true_x)-1, max(true_x)+1, 100);
  plot(x, f_transed(x));
  plot(eval_x, approxi_y, '^');
  legend('Given Data', 'Estimated curve', 'x=2.6', 'location', 'best');
  ```
### Result
  ```
  >> HW06_Q4
  a=2.172916834412650
  b=0.499435719499467
  y =
    function_handle with value:
      @(x)a*log(x)+b

  x =
    function_handle with value:
      @(y)log((y-b)/a)

  approxi_y = 2.575682623873542
  ```
  |  |  |
  | -------- | -------- |
  | ![](https://i.imgur.com/0eJ7odn.png)     | ![](https://i.imgur.com/fRaCbX3.png)     |

---
## Q5 Least-squares Regression with Matrix
The following data provided:
$\ 　　\begin{array} {|r|r|}\hline x & 1 & 2 & 3 & 3 & 4 \\ \hline y & 2.2 & 2.8 & 3.6 & 4.5 & 5.5 \\ \hline  \end{array}$
You have to use the least-squares regression to fit this data with the following model:
  $\ 　　y=a+bx+\dfrac{c}{x}$
Determine the coefficients by setting up and solving:
  $\ 　　\left( Z^TZ\right)A=Z^TY$
### Calculation
- Declartions
  We have given data
  $\ 　X=\left[ 1　2　3　4　5\right]^T$
  $\ 　Y=\left[ 2.2　2.8　3.6　4.5　5.5\right]^T$
  Define $Z$ as:
  $\ 　Z=\begin{bmatrix}
    1&X_1&\frac{1}{X_1}\\ \:1&X_2&\frac{1}{X_2}\\ \:1&X_3&\frac{1}{X_3}\\ \:1&X_4&\frac{1}{X_4}\\ \:1&X_5&\frac{1}{X_5}
    \end{bmatrix}$
  ```c
  true_x = [1 2 3 4 5]';
  true_y = [2.2 2.8 3.6 4.5 5.5]';
  f = @(x) a + b.*x + c./x;
  
  n = length(true_x);
  Z = [ones(n,1) true_x 1./true_x];
  ```
- Solve
  I define the coefficients $a, b, c$ as $A$
  $\ 　A=[a　b　c]^T$
  Solve for $A$:
  $\ 　A=\left( Z^TZ\right)^{-1}Z^TY$
  ```c
  A = inv(Z'*Z)*Z'*true_y;  % //(Z' * Z) * A = Z' * true_y, solve for A

  a = A(1)
  b = A(2)
  c = A(3)
  f = @(x) a + b.*x + c./x;
  ```
- Verify Coefficient of Determinat
  We have solved the matrix $A$. The coefficients $a, b, c$ are known now.
  $\ 　　f(x)=y=a+bx+\dfrac{c}{x}$
  $\ 　　S_t=\sum{\left(Y_i-Y_{avg}\right)^2}$
  $\ 　　S_r=\sum{\left(Y_i-f(X_i)\right)^2}$
  $\ 　　　　\text{where}\ X_i\in X, Y_i\in Y\ \text{and}\ Y_{avg}=\dfrac{\sum{Y_i}}{n}$
  $\ 　　r^2=\dfrac{S_t-S_r}{S_t}$
  ```c
  y_avg = mean(true_y);
  S_t = sum( (true_y-y_avg).^2 );
  S_r = sum( (true_y-f(true_x)).^2 );
  r_square = (S_t-S_r) / S_t
  ```
### Code
  ```c=
  clear all
  close all
  format long

  % //Given data, regresion 2 to y = a + b*x + c/x, find a, b and c
  true_x = [1 2 3 4 5]';
  true_y = [2.2 2.8 3.6 4.5 5.5]';

  % //Least-squares Regression by matrix
  n = length(true_x);
  Z = [ones(n,1) true_x 1./true_x];
  A = inv(Z'*Z)*Z'*true_y;  % //(Z' * Z) * A = Z' * true_y, solve for A

  % //y = a + b*x + c/x
  a = A(1)
  b = A(2)
  c = A(3)
  f = @(x) a + b.*x + c./x;

  % //Verify Coefficient of Determinat
  y_avg = mean(true_y);
  S_t = sum( (true_y-y_avg).^2 );
  S_r = sum( (true_y-f(true_x)).^2 );
  r_square = (S_t-S_r) / S_t

  hold on;
  plot(true_x, true_y, 'o');
  x = linspace(min(true_x)-1, max(true_x)+1, 100);
  plot(x, f(x));
  legend('Given Data', 'Estimated curve');
  title('Least-squares Regression');
  ```
### Result
  ```
  >> HW06_Q5
  a =
     0.374496644295302

  b =
     0.986442953020126

  c =
     0.845637583892596

  r_square =
     0.999602029264372     
  ```
  ![](https://i.imgur.com/ThoYiGp.png =500x)

  
  