---
tags: homeworks, Numerical Methods
---
# Numerical Methods Homework-7
## B10702026 林琨霖
[Online link of this homework](https://hackmd.io/khMidJcgSiqgyQNwAe5Bag?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 Fit $f(t)=\text{sin} ^2(t)$
### (a) Not-a-knot
Cubic spline with not-a-knot end conditions.
- Code
```c=
clear
clear all
close all
format long
%// Ref: http://www.cs.tau.ac.il/~turkel/notes/numeng/spline_note.pdf
%// Given data --------------------------------------------------------------
f = @(t) sin(t).^2;
t_start = 0;
t_end = 2*pi;
x = linspace(t_start, t_end, 8);
y = f(x);
%// Evaluate spacing h and slopes d -----------------------------------------
%// lenght of x and y must match
n = length(x); %// Given points counts
h = zeros(1,n-1); %// x-spacing between points
d = zeros(1,n-1); %// delta_y/delta_x
for i = 1:(n-1)
h(i) = x(i+1)-x(i);
d(i) = (y(i+1)-y(i)) / h(i);
end
%// Not-a-knot Spline -------------------------------------------------------
A = zeros(n,n);
b = zeros(n,1);
A(1,:) = [h(2) -(h(1)+h(2)) h(1) linspace(0,0,n-3)];
A(n,:) = [linspace(0,0,n-3) h(n-1) -(h(n-2)+h(n-1)) h(n-2)];
b(1) = 0; b(n) = 0;
for i = 2:(n-1)
a = linspace(0,0,n);
a(i-1) = h(i-1);
a(i) = 2*(h(i-1)+h(i));
a(i+1) = h(i);
A(i,:) = a;
b(i) = 6*( d(i)-d(i-1) );
end
m = inv(A)*b; %// A*m = b
%// Build peice wise splines ------------------------------------------------
s=zeros(n-1,4);
s(:,1) = y(1:1:n-1)'; %// yk
s(:,2) = d' - h'./6.*( 2.*m(1:1:n-1) + m(2:1:n) ); %// dk-hk / 6*(2mk+m(k+1))
s(:,3) = m(1:1:n-1)./2; %// mk/2
s(:,4) = (m(2:1:n)-m(1:1:n-1)) ./ (6.*h'); %// (m(k+1)-mk) / (6hk)
s_func = @(i,x_) s(i,1) + s(i,2).*(x_-x(i)) + s(i,3).*(x_-x(i)).^2 + s(i,4).*(x_-x(i)).^3;
%// Plot --------------------------------------------------------------------
x_plots = linspace(x(1), x(n), 200);
%// The first point is skipped in the following for
y_plots = s_func(1, x(1));
for i = 1:(n-1)
y_plots = [y_plots s_func(i, x_plots(x(i)< x_plots & x_plots<=x(i+1)) )];
end
hold on;
plot(x, y, 'og');
plot(x_plots, y_plots);
plot(x_plots, f(x_plots));
legend("Eight points", "Spline", 'f(t)', 'location', 'best');
title("Not-a-not Spline with original function");
xlabel("x"); ylabel("y");
```
- Result
```
>> HW07_Q1_a
>> s
s =
0 0.224904124663937 0.846599063847394 -0.377088971299779
0.611260466978157 0.833273829513334 -0.168823743568421 -0.377088971299779
0.950484433951210 -0.381239226817294 -1.184246550984235 0.738540126662379
0.188255099070633 -0.722108222638808 0.804489651561884 0.000000000000000
0.188255099070633 0.722108222638808 0.804489651561884 -0.738540126662379
0.950484433951210 0.381239226817295 -1.184246550984235 0.377088971299778
0.611260466978157 -0.833273829513334 -0.168823743568421 0.377088971299778
>> m
m =
1.693198127694788
-0.337647487136841
-2.368493101968471
1.608979303123768
1.608979303123769
-2.368493101968470
-0.337647487136842
1.693198127694786
```
| | |
| ------------------------------------ | ------------------------------------ |
|  |  |
### (b) Derivative end conditions
Cubic spline with with derivative end conditions equal to the exact values calculated
with differentiation.
The end condition described above matches to the end condition of Clamped Spline.
- Code
```c=
clear
clear all
close all
format long
%// Ref: http://www.cs.tau.ac.il/~turkel/notes/numeng/spline_note.pdf
%// Given data --------------------------------------------------------------
f = @(t) sin(t).^2;
f_diff = @(t) 2.*sin(t).*cos(t);
t_start = 0;
t_end = 2*pi;
x = linspace(t_start, t_end, 8);
y = f(x);
%// Evaluate spacing h and slopes d -----------------------------------------
%// lenght of x and y must match
n = length(x); %// Given points counts
h = zeros(1,n-1); %// x-spacing between points
d = zeros(1,n-1); %// delta_y/delta_x
for i = 1:(n-1)
h(i) = x(i+1)-x(i);
d(i) = (y(i+1)-y(i)) / h(i);
end
%// Clamped Spline -------------------------------------------------------
A = zeros(n,n);
b = zeros(n,1);
A(1, 1) = 2 * h(1);
A(1, 2) = h(1);
A(n, n-1) = h(n - 1);
A(n, n) = 2 * h(n - 1);
b(1) = 6 * (d(1) - f_diff(x(1)));
b(n) = 6 * (f_diff(x(n-1)) - d(n-1));
for i = 2:(n-1)
a = linspace(0,0,n);
a(i-1) = h(i-1);
a(i) = 2*(h(i-1)+h(i));
a(i+1) = h(i);
A(i,:) = a;
b(i) = 6*( d(i)-d(i-1) );
end
m = inv(A)*b; %// A*m = b
%// Build peice wise splines ------------------------------------------------
s=zeros(n-1,4);
s(:,1) = y(1:1:n-1)'; %// yk
s(:,2) = d' - h'./6.*( 2.*m(1:1:n-1) + m(2:1:n) ); %// dk-hk / 6*(2mk+m(k+1))
s(:,3) = m(1:1:n-1)./2; %// mk/2
s(:,4) = (m(2:1:n)-m(1:1:n-1)) ./ (6.*h'); %// (m(k+1)-mk) / (6hk)
s_func = @(i,x_) s(i,1) + s(i,2).*(x_-x(i)) + s(i,3).*(x_-x(i)).^2 + s(i,4).*(x_-x(i)).^3;
%// Plot --------------------------------------------------------------------
x_plots = linspace(x(1), x(n), 200);
%// The first point is skipped in the following for
y_plots = s_func(1, x(1));
for i = 1:(n-1)
y_plots = [y_plots s_func(i, x_plots(x(i)< x_plots & x_plots<=x(i+1)) )];
end
hold on;
plot(x, y, 'og');
plot(x_plots, y_plots);
plot(x_plots, f(x_plots));
legend("Eight points", "Spline", 'f(t)', 'location', 'best');
title("Clamped Spline with original function");
xlabel("x"); ylabel("y");
```
- Result
> 
> The result is quite similar to Q1(a)
### \(c\) PCHIP
Piecewise cubic Hermite interpolation polymonial (PCHIP).
- Code
```c=
clear
clear all
close all
format long
%// Ref: http://www.cs.tau.ac.il/~turkel/notes/numeng/spline_note.pdf
%// Given data --------------------------------------------------------------
f = @(t) sin(t).^2;
f_diff = @(t) 2.*sin(t).*cos(t);
t_start = 0;
t_end = 2*pi;
x = linspace(t_start, t_end, 8);
y = f(x);
n = length(x);
%// Plot --------------------------------------------------------------------
x_plots = linspace(x(1), x(n), 200);
y_plots = interp1(x, y, x_plots, 'pchip');
hold on;
plot(x, y, 'og');
plot(x_plots, y_plots);
plot(x_plots, f(x_plots));
legend("given points", "PCHIP", 'f(t)', 'location', 'best');
title("PCHIP with original function");
xlabel("x"); ylabel("y");
```
- Result
## Q2 Newton’s interpolation
Calculate $f(4)$ using Newton's interpolating polynomials of order 1 through 4. Choose your base points to attain good accuracy. What do your resultsindicate regarding the order of the polynomial used to generate the data in the table?
| $x$ | 1 | 2 | 3 | 5 | 7 | 8 |
|:----:|:-:|---|----|----|-----|:---:|
|$f(x)$| 3 | 6 | 19 | 99 | 291 | 444 |
- Code
```c=
clear
clear all
close all
format long
%// Given data
x = [1 2 3 5 7 8];
y = [3 6 19 99 291 444];
x_tar = 4; %// Estimate f(x=4)
n = length(x);
x_plots = linspace(x(1), x(n), 200);
y_plots = zeros(n-1,200);
highest_order = 4;
%// 1st~5th order Newton interpolations
for order = 1:highest_order
%// Select base points near x_tar
switch order
case 1, index = 3:4;
case 2, index = 2:4;
case 3, index = 2:5;
case 4, index = 1:5;
otherwise, index = 1:n;
end
y_plots(order,:) = Newtint(x(index), y(index), x_plots);
y_eval = Newtint(x(index), y(index), x_tar);
fprintf("No.%d order with base points:\n", order);
for i=1:order, fprintf("(%d,%d) ", x(i), y(i)); end
fprintf("\nf(%d) = %.15f\n\n", x_tar, y_eval);
end
%// Plot ----------------------------------------
hold on;
plot(x, y, 'xb'); %// Given data
for i = 1:highest_order
plot(x_plots, y_plots(i,:));
end
legend('Given data','1_{st} order','2_{nd} order','3_{rd} order',...
'4_{th} order', '5_{th} order', ...
'location', 'best');
title("Newton interpolation of different orders");
xlabel("x"); ylabel("y");
```
> The function Newtint() is in the appendix.
- Result
```
>> HW07_Q2
1st order with base points:
(3,19) (5,99)
f(4) = 59.0000000
2nd order with base points:
(2,6) (3,19) (5,99)
f(4) = 50.0000000
3rd order with base points:
(2,6) (3,19) (5,99) (7,291)
f(4) = 48.0000000
4th order with base points:
(1,3) (2,6) (3,19) (5,99) (7,291)
f(4) = 48.0000000
```
| ith order | $f_i(4)$ |
|:---------:|:--------:|
| 1 | 59 |
| 2 | 50 |
| 3 | 48 |
| 4 | 48 |
> 
> Note that the curve of the 3rd order is overlapped with the 4th order.
- Discussion
Define $f_i(x)$ is the interpolated function, where $i$ is the order of Newton's interpolating polynomials.
Since the original function $f(x)$ is not given, I can say $f_4(x)$ yields the closest value to $f(x)$. So $f_4(4)=48$ is the most accurate value among $f_3(4)$, $f_2(4)$ and $f_1(4)$.
- Base points
For $f_1(x)$, I choose the base points $(3,19)$ and $(5,99)$. Because $x$ of these two points are the cloesest to $4$, as the required approximation of $f(4)$. Similarly, for $f_2(x)$, three points near by $x=4$ are choosen: $(2,6)$, $(3,19)$ and $(5,99)$ are choosen.
- Marginal utility
Recall the graph result, the curve of the 3rd order is overlapped with the 4th order. I wonder if the curve of the 5th order overlapped. Since there are 6 points given, so the highest order can be calculated is 5. Here is the result:
```
5th order with base points:
(1,3) (2,6) (3,19) (5,99) (7,291) (8,444)
f(4) = 48.0000000
```
It turns out that the curve of 5th order is overlapped too. And the interpolated function yields same result $f_5(4)=48$.
As marginal utility described in Economics, with each increasing order, the improvement of the result is smaller and smaller.
---
## Q3 Lagrange interpolation
Repeat Q2 by using Lagrange polynomials of order 1 through 3.
- Code
```c=
clear
clear all
close all
format long
%// Given data
x = [1 2 3 5 7 8];
y = [3 6 19 99 291 444];
x_tar = 4; %// Estimate f(x=4)
n = length(x);
x_plots = linspace(x(1), x(n), 200);
y_plots = zeros(n-1,200);
%// 1st~3th order Lagrange interpolations
for order = 1:3
%// Select base points near x_tar
switch order
case 1, index = 3:4;
case 2, index = 2:4;
case 3, index = 2:5;
case 4, index = 1:5;
otherwise, index = 1:n;
end
y_plots(order,:) = Lagrange(x(index), y(index), x_plots);
y_eval = Lagrange(x(index), y(index), x_tar);
fprintf("No.%d order with base points:\n", order);
for i=1:order, fprintf("(%d,%d) ", x(i), y(i)); end
fprintf("\nf(%d) = %.15f\n\n", x_tar, y_eval);
end
%// Plot ----------------------------------------
hold on;
plot(x, y, 'xb'); %// Given data
for i = 1:3
plot(x_plots, y_plots(i,:));
end
legend('Given data','1_{st} order','2_{nd} order','3_{rd} order',...
'location', 'best');
title("Lagrange interpolation of different orders");
xlabel("x"); ylabel("y");
```
> The function Lagrange() is in the appendix.
- Result
> The result is quite similar to Q2
> 
---
## Q4 Runge's function
$f(x)=\dfrac{1}{1+25x^2}$
### (a) Plot
- Code
```c=
f = @(x) 1 ./ (1+25.*x.^2);
x_plot = linspace(-1, 1, 100);
y_plot = f(x_plot);
hold on;
plot(x_plot, y_plot);
leg1 = legend('$\frac{1}{1+25x^2}$','Interpreter','latex');
set(leg1,'FontSize',17);
title("Runge's function");
```
- Result
### (b) Lagrange interpolating
Generate and plot e forth-order Lagrange interpolating polynomial using equispaced
function values corresponding to $x =-1,-0.5,0,0.5,\text{ and } 1$.
- Code
```c=
clear
clear all
close all
format long
%// Given data
f = @(x) 1 ./ (1+25.*x.^2);
x = [-1 -0.5 0 0.5 1];
y = f(x);
%// 4th order of Lagrange
n = length(x);
x_plots = linspace(x(1), x(n), 200);
y_plots = Lagrange(x, y, x_plots);
%// Plot ----------------------------------------
hold on;
plot(x, y, 'xb'); %// Given data
plot(x_plots, y_plots);
plot(x_plots, f(x_plots));
legend('Given data','4_{th} order','Runge’s function','location', 'best');
title("Lagrange of 4th order fitting Runge’s function");
xlabel("x"); ylabel("y");
```
```c=
function yint = Lagrange(x,y,xx)
n = length(x);
if length(y)~=n, error('x and y must be same length'); end
s = 0;
for i = 1:n
product = y(i);
for j = 1:n
if i ~= j
product = product .* (xx-x(j)) ./ (x(i)-x(j));
end
end
s = s+product;
end
yint = s;
```
- Result
### \(c\) Continued from (b)
Use five points from (b) to estimate $f(0.8)$ with 1st through 4th order Newton interpolating polynomials.
- Code
```c=
clear all
close all
format long
%// Given data
f = @(x) 1 ./ (1+25.*x.^2);
x = [-1 -0.5 0 0.5 1];
y = f(x);
x_tar = 0.8; %// Estimate f(x=0.8)
n = length(x);
x_plots = linspace(x(1), x(n), 200);
y_plots = zeros(n-1,200);
%// 1st~4th order Newton interpolations
for order = 1:4
%// Select base points near x_tar
switch order
case 1, index = 4:5;
case 2, index = 3:5;
case 3, index = 2:5;
case 4, index = 1:5;
otherwise, index = 1:n;
end
y_plots(order,:) = Newtint(x(index), y(index), x_plots);
y_eval = Newtint(x(index), y(index), x_tar);
fprintf("No.%d order with base points:\n", order);
for i=1:order, fprintf("(%d,%d) ", x(i), y(i)); end
fprintf("\nf(%d) = %.15f\n\n", x_tar, y_eval);
end
%// Plot ----------------------------------------
hold on;
plot(x, y, '^b'); %// Given data
for i = 1:4
plot(x_plots, y_plots(i,:));
end
plot(x_plots, f(x_plots));
legend('Given data','1_{st} order','2_{nd} order','3_{rd} order',...
'4_{th} order','Runge’s function', 'location', 'best');
title("Newton interpolation of different orders");
xlabel("x"); ylabel("y");
```
- Result
```
>> HW07_Q4_c
1st order with base points:
(5.000000e-01,1.379310e-01) (1,3.846154e-02)
f(8.000000e-01) = 0.0782493
2nd order with base points:
(0,1) (5.000000e-01,1.379310e-01) (1,3.846154e-02)
f(8.000000e-01) = -0.0132626
3rd order with base points:
(-5.000000e-01,1.379310e-01) (0,1) (5.000000e-01,1.379310e-01) (1,3.846154e-02)
f(8.000000e-01) = -0.1724138
4th order with base points:
(-1,3.846154e-02) (-5.000000e-01,1.379310e-01) (0,1) (5.000000e-01,1.379310e-01) (1,3.846154e-02)
f(8.000000e-01) = -0.3793103
```
| ith order | $f_i(0.8)$ | ｜$f(x)-f_i(x)$｜ |
|:---------:|:------------------:| ----------------- |
| original | 0.058823529411765 | 0 |
| 1 | 0.078249336870027 | 0.019425807458262 |
| 2 | -0.013262599469496 | 0.072086128881261 |
| 3 | -0.172413793103449 | 0.231237322515214 |
| 4 | -0.379310344827586 | 0.438133874239351 |
> 
> 
> 
> From this graph and the table above, the error($\big|f(x)-f_i(x)\big|$) is greater with higher order around $x=0.8$.
> But when viewing at the range of $-1\leq x\leq 1$ (first picture of this section), there is not a curve that fits Runge's function well within the range.
### (d) Cubic spline
Generate and plot a cubic spline using the five points from (b).
- Code
```c=
clear all
close all
format long
%// Given data
f = @(x) 1 ./ (1+25.*x.^2);
x = [-1 -0.5 0 0.5 1];
y = f(x);
x_tar = 0.8; %// Estimate f(x=0.8)
n = length(x);
x_plots = linspace(x(1), x(n), 200);
%// Splines -------------------------------
%// cnd=1: Not a knot; cnd=2: Natural
for cnd=1:2
h = zeros(1,n-1); %// x-spacing between points
d = zeros(1,n-1); %// delta_y/delta_x
for i = 1:(n-1)
h(i) = x(i+1)-x(i);
d(i) = (y(i+1)-y(i)) / h(i);
end
%// Build and solve for A*m=b
A = zeros(n,n);
b = zeros(n,1);
switch cnd
case 1 %// Not-a-knot Spline end conditions
A(1,:) = [h(2) -(h(1)+h(2)) h(1) linspace(0,0,n-3)];
A(n,:) = [linspace(0,0,n-3) h(n-1) -(h(n-2)+h(n-1)) h(n-2)];
case 2 %// Natural Spline end conditions
A(1,:) = [1 linspace(0,0,n-1)];
A(n,:) = [linspace(0,0,n-1) 1];
end
b(1) = 0; b(n) = 0;
for i = 2:(n-1)
a = linspace(0,0,n);
a(i-1) = h(i-1);
a(i) = 2*(h(i-1)+h(i));
a(i+1) = h(i);
A(i,:) = a;
b(i) = 6*( d(i)-d(i-1) );
end
m = A\b;
%// Build peice wise splines
s = zeros(n-1,4);
s(:,1) = y(1:1:n-1)'; %// yk
s(:,2) = d' - h'./6.*( 2.*m(1:1:n-1) + m(2:1:n) ); %// dk-hk / 6*(2mk+m(k+1))
s(:,3) = m(1:1:n-1)./2; %// mk/2
s(:,4) = (m(2:1:n)-m(1:1:n-1)) ./ (6.*h'); %// (m(k+1)-mk) / (6hk)
s_func = @(i,x_) s(i,1) + s(i,2).*(x_-x(i)) + s(i,3).*(x_-x(i)).^2 + s(i,4).*(x_-x(i)).^3;
%// Fit the points
y_pts = s_func(1, x(1));
%// Evaluates y at x1~x2, x2~x3, x3~4...
for i = 1:(n-1)
y_pts = [y_pts s_func(i, x_plots(x(i) < x_plots & x_plots <= x(i+1)) )];
end
y_plots(cnd,:) = y_pts;
end
%// Plot --------------------------------------------------------------------
x_plots = linspace(x(1), x(n), 200);
hold on;
plot(x, y, '^b'); %// Given data
plot(x_plots, y_plots(1,:)); %// Not a knot Splines
plot(x_plots, y_plots(2,:)); %// Natural Splines
plot(x_plots, f(x_plots)); %// Original function
legend("Given points", "Not-a-knot Splines","Natural Splines","Runge’s function", 'location', 'best');
title("Splines fitting Runge’s function");
xlabel("x"); ylabel("y");
```
- Result
> 
> Recall the result of Q4(b), the shape of esitimated curve is quite similar.
> Both of the splines result osscilate over Runge's function. Yet, Natrual splines seems to have smaller error within the range $-1\leq x\leq 1$.
### (e) Discusssion
I overlapped the result of Q4(b\) and Q4(e):
It's obvious that Natural Splines yields a better result than Lagrange interpolation. But the error of Natural Splines is still not neglectable in this scope. I think it's because Runge's function has massive changes only around $x=0$, while cubic splines and Lagrange interpolation seem to suppose that
> $y$ has greater changes between the given points which has minor changes in $y$.
For example, the changes in $y$ between the first given point($x=-1$) and the second one($x=-0.5$) is small. It results in big swings in $y$ in this range($-1\leq x\leq -0.5$) in both Natural Splines and Lagrange interpolation.
What's more, the Spline method suppose the end condition as the terminal points ($x=-1,1$ in this case) are the inflection points.
$$m_1=0,\ m_n=0$$
where $n$ is the counts of given points. This is not the case in Runge's function.
I also tried different *given points*(or *sampling points*) to fit the Runge's function, but none of the results matches the original function well.
> Trying differen given points
> |  |  |
> |:------------------------------------:|:-------------------------------------------------------------------:|
> |  |  |
> |  | Fitting with 9 given points<br> |
Natural splines only fits the Runge's funtion well when the counts of given points is great enough, while Lagrange swings more severely. But increasing the *given points* is not an available option in the most of applications, not to mention I know what the target function is before genrating a approximation.
## Appendix
```c=
function yint = Newtint(x,y,xx)
%// Newtint:Newton interpolating polynomial
%// yint=Newtint(x,y,xx):Uses an (n-1)-order Newton
%// interpolating polynomial based on n data points (x,y)
%// to determine a value of the dependent variable (yint)
%// Aat a given value of the independent variable,xx.
%// input:
%// x=independent variable
%// y=dependent variable
%// xx=value of independent variable at which interpolation is calculated
%// output:
%// yint=interpolated value of dependent variable
%// compute the finite divided differences in the form of a difference table
n=length(x);
if length(y)~=n, error('x and y must be same length');
end
b=zeros(n,n);
%// assign dependent variables to the first columm of b.
b(:,1)=y(:); %the (:)ensures that y is a column vector.
for j=2:n
for i=1:n-j+1
b(i,j)=(b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i));
end
end
%// use the finite divided differences to interpolate
xt=1;
yint=b(1,1);
for j=1:n-1
xt=xt.*(xx-x(j));
yint=yint+b(1,j+1).*xt;
end
```
```c=
function yint = Lagrange(x,y,xx)
%// Lagrange interpolation. Uses an (n - 1)-order Lagrange
%// interpolating polynomial based on n data points (x, y)
%// to determine a value of the dependent variable (yint)
%// at a given value of the independent variable, xx.
%// input:
%// x = independent variable
%// y = dependent variable
%// xx = value of independent variable at which the
%// interpolation is calculated
%// output:
%// yint = interpolated value of dependent variable
n = length(x);
if length(y)~=n, error('x and y must be same length'); end
s = 0;
for i = 1:n
product = y(i);
for j = 1:n
if i ~= j
product = product .* (xx-x(j)) ./ (x(i)-x(j));
end
end
s = s+product;
end
yint = s;
```
Sign in with Wallet
Connect another wallet