---
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.
> 
## 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() | x-y-f() |
| x-f()<br> | y-f()<br> |
| x-y (Origin at bottom-right corner)<br> | x-y (Origin at upper-left corner)<br> |
---
## 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> | x-y<br> |
| --------------------------------------------------------- | --------------------------------------------------------- |
| x-f(x,y)<br> | y-f(x,y)<br> |
|  |  |
|  |  |
---
## 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.
### 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 |  | Two lines overlapped  |
| Location<br>v.s.<br>Time |  |  |
---
## 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
```
| | |
| -------- | -------- |
|  |  |
---
## 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
```
Sign in with Wallet
Connect another wallet