notion：
https://www.notion.so/hw1-Numerical-Analysis-bd1529dfaa7d4c6c8216da1e4fbbc684
In this homework we will use Lagrange polynomial interpolation and cubic spline interpolation to interpolate (N + 1 = 11) data points (xi , yi), where i = 0, 1, 2, ..., 10. The independent variable x is in the range of [−1, 1]. There are two sets of (N + 1 = 11) data points.
Part A. (25%) Please refer to the file hw1AB.dat for the 11 data points. Here xi , where i = 0, 1, 2, ..., 10, are uniformly distributed in [−1, 1]. A.1 Plot the Lagrange polynomial Lj (x), where j = 0, 1, 2, ..., 10. Note that Lj (xi) = 0 when i 6= j and Lj (xi) = 1 when i = j. A.2 Plot the interpolating polynomial that goes through the 11 data points P(x) = X 10 j=0 yjLj (x)
Part B. (25%) Please refer to the file hw1AB.dat for the 11 data points. Here xi , where i = 0, 1, 2, ..., 10, are uniformly distributed in [−1, 1]. B.1 Here we use cubic spline interpolation and we should assume the second derivative at the points, g ′′(xi) where i = 0, 1, 2, ..., 10, as unknowns. Let’s use free run-out condition for g ′′(x0) = g ′′(x10) = 0. What are the values of g ′′(xi) where i = 1, 2, ..., 9? B.2 Plot the cubic spline interpolation for the whole range of x ∈ [−1, 1].
Part C. (25%) Please refer to the file hw1CD.dat for the 11 data points. Here xi , where i = 0, 1, 2, ..., 10, are non-uniformly distributed in [−1, 1]. C.1 Plot the Lagrange polynomial Lj (x), where j = 0, 1, 2, ..., 10. Note that Lj (xi) = 0 when i 6= j and Lj (xi) = 1 when i = j. C.2 Plot the interpolating polynomial that goes through the 11 data points P(x) = X 10 j=0 yjLj (x) 1
Part D. (25%) Please refer to the file hw1CD.dat for the 11 data points. Here xi , where i = 0, 1, 2, ..., 10, are non-uniformly distributed in [−1, 1]. D.1 Here we use cubic spline interpolation and we should assume the second derivative at the points, g ′′(xi) where i = 0, 1, 2, ..., 10, as unknowns. Let’s use free run-out condition for g ′′(x0) = g ′′(x10) = 0. What are the values of g ′′(xi) where i = 1, 2, ..., 9? D.2 Plot the cubic spline interpolation for the whole range of x ∈ [−1, 1].
hw1AB.dat
xj yj
-1.0000 0.0385
-0.8000 0.0588
-0.6000 0.1000
-0.4000 0.2000
-0.2000 0.5000
0.0000 1.0000
0.2000 0.5000
0.4000 0.2000
0.6000 0.1000
0.8000 0.0588
1.0000 0.0385
hw2CD.dat
xj yj
-1.0000 0.0385
-0.9511 0.0424
-0.8090 0.0576
-0.5878 0.1038
-0.3090 0.2952
0.0000 1.0000
0.3090 0.2952
0.5878 0.1038
0.8090 0.0576
0.9511 0.0424
1.0000 0.0385
filename = "hw1AB.dat";
[datax, datay] = textread(filename, "%f %f", 'headerlines', 1);
x0 = datax;
y0 = datay;
x = linspace(-1, 1, 500)
d = size(x0, 1);
y = zeros(d, 500);
y = y+1;
for j = 1:d
for i = 1:d
if j == i
continue;
#skip and continue the loop
endif
y(j, :) .*= ((x-x0(i))/(x0(j)-x0(i)));
endfor
endfor
for k = 1:1:d
subplot(6, 2, k)
plot(x, y(k, :))
set(gca,'FontSize',10);
xlabel("X",'FontSize',15);
ylabel("Lj(x)",'FontSize',15);
title(x0(k),'FontSize',15);
grid on;
end
* the function should be written in another file in the same folder
B.2
filename = 'hw1AB.dat'
[xj,yj]=textread(filename,'%f %f','HeaderLines',1)
n=length(xj);
g=zeros(n-2,n-2);
f=zeros(n-2,1);
G=zeros(n-2,1);
h=zeros(n-1,1);
for o=1:n-1;
h(o)=xj(o+1)-xj(o);
endfor;
for i=1:n-2;
for j=1:n-2;
if j==i-1;
g(i,j)=h(i)./6;
endif;
if j==i;
g(i,j)=(h(i+1).+h(i))./3;
endif;
if j==i+1;
g(i,j)=h(i+1)./6;
endif;
endfor;
endfor;
for k=1:n-2;
f(k,1)=((yj(k+2).-yj(k+1))./h(k+1)).-((yj(k+1).-yj(k))./h(k));
endfor;
G=g\f;
x=-1:0.0001:-0.8001;
x1=-0.8:0.0001:-0.6001;
x2=-0.6:0.0001:-0.4001;
x3=-0.4:0.0001:-0.2001;
x4=-0.2:0.0001:-0.0001;
x5=0:0.0001:0.1999;
x6=0.2:0.0001:0.3999;
x7=0.4:0.0001:0.5999;
x8=0.6:0.0001:0.8;
x9=0.8001:0.0001:1;
c=1;
c1=1;
c2=1;
c3=1;
c4=1;
c5=1;
c6=1;
c7=1;
c8=1;
c9=1;
c .*=G(1)./6.*(((x.-xj(1)).^3)./h(1).-(h(1).*(x.-xj(1)))) .+ ((yj(1).*(xj(2).-x).+yj(2).*(x.-xj(1)))./h(1));
c1 .*=(G(1)./6).*(((xj(3).-x1).^3)./h(2).-(h(2).*(xj(3).-x1))).+(G(2)./6).*(((x1.-xj(2)).^3)./h(2).-(h(2).*(x1.-xj(2)))).+((yj(2).*(xj(3).-x1).+yj(3).*(x1.-xj(2)))./h(2));
c2 .*=(G(2)./6).*(((xj(4).-x2).^3)./h(3).-(h(3).*(xj(4).-x2))).+(G(3)./6).*(((x2.-xj(3)).^3)./h(3).-(h(3).*(x2.-xj(3)))).+((yj(3).*(xj(4).-x2).+yj(4).*(x2.-xj(3)))./h(3));
c3 .*=(G(3)./6).*(((xj(5).-x3).^3)./h(4).-(h(4).*(xj(5).-x3))).+(G(4)./6).*(((x3.-xj(4)).^3)./h(4).-(h(4).*(x3.-xj(4)))).+((yj(4).*(xj(5).-x3).+yj(5).*(x3.-xj(4)))./h(4));
c4 .*=(G(4)./6).*(((xj(6).-x4).^3)./h(5).-(h(5).*(xj(6).-x4))).+(G(5)./6).*(((x4.-xj(5)).^3)./h(5).-(h(5).*(x4.-xj(5)))).+((yj(5).*(xj(6).-x4).+yj(6).*(x4.-xj(5)))./h(5));
c5 .*=(G(5)./6).*(((xj(7).-x5).^3)./h(6).-(h(6).*(xj(7).-x5))).+(G(6)./6).*(((x5.-xj(6)).^3)./h(6).-(h(6).*(x5.-xj(6)))).+((yj(6).*(xj(7).-x5).+yj(7).*(x5.-xj(6)))./h(6));
c6 .*=(G(6)./6).*(((xj(8).-x6).^3)./h(7).-(h(7).*(xj(8).-x6))).+(G(7)./6).*(((x6.-xj(7)).^3)./h(7).-(h(7).*(x6.-xj(7)))).+((yj(7).*(xj(8).-x6).+yj(8).*(x6.-xj(7)))./h(7));
c7 .*=(G(7)./6).*(((xj(9).-x7).^3)./h(8).-(h(8).*(xj(9).-x7))).+(G(8)./6).*(((x7.-xj(8)).^3)./h(8).-(h(8).*(x7.-xj(8)))).+((yj(8).*(xj(9).-x7).+yj(9).*(x7.-xj(8)))./h(8));
c8 .*=(G(8)./6).*(((xj(10).-x8).^3)./h(9).-(h(9).*(xj(10).-x8))).+(G(9)./6).*(((x8.-xj(9)).^3)./h(9).-(h(9).*(x8.-xj(9)))).+((yj(9).*(xj(10).-x8).+yj(10).*(x8.-xj(9)))./h(9));
c9 .*=(G(9)./6).*(((xj(11).-x9).^3)./h(10).-(h(10).*(xj(11).-x9))) .+ ((yj(10).*(xj(11).-x9).+yj(11).*(x9.-xj(10)))./h(10));
y=c;
y1=c1;
y2=c2;
y3=c3;
y4=c4;
y5=c5;
y6=c6;
y7=c7;
y8=c8;
y9=c9;
plot(x,y,x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,x7,y7,x8,y8,x9,y9);
set(gca,'FontSize',20);
xlabel("X",'FontSize',20);
ylabel("g(x)",'FontSize',20);
title("g(x) with 11 datas in hw1AB",'FontSize',20);
grid on;
* the function should be written in another file in the same folder
filename = 'hw1CD.dat'
[xj,yj]=textread(filename,'%f %f','HeaderLines',1)
n=length(xj);
g=zeros(n-2,n-2);
f=zeros(n-2,1);
G=zeros(n-2,1);
h=zeros(n-1,1);
x=-1:0.0001:-0.9512;
x1=-0.9511:0.0001:-0.8091;
x2=-0.8090:0.0001:-0.5879;
x3=-0.5878:0.0001:-0.3091;
x4=-0.3090:0.0001:-0.0001;
x5=0:0.0001:0.3089;
x6=0.3090:0.0001:0.5877;
x7=0.5878:0.0001:0.8089;
x8=0.8090:0.0001:0.9511;
x9=0.9512:0.0001:1;
c=1;
c1=1;
c2=1;
c3=1;
c4=1;
c5=1;
c6=1;
c7=1;
c8=1;
c9=1;
for o=1:n-1;
h(o)=xj(o+1)-xj(o);
endfor;
for i=1:n-2;
for j=1:n-2;
if j==i-1;
g(i,j)=h(i)./6;
endif;
if j==i;
g(i,j)=(h(i+1).+h(i))./3;
endif;
if j==i+1;
g(i,j)=h(i+1)./6;
endif;
endfor;
endfor;
for k=1:n-2;
f(k,1)=((yj(k+2).-yj(k+1))./h(k+1)).-((yj(k+1).-yj(k))./h(k));
endfor
G=g\f;
c .*=G(1)./6.*(((x.-xj(1)).^3)./h(1).-(h(1).*(x.-xj(1)))) .+ ((yj(1).*(xj(2).-x).+yj(2).*(x.-xj(1)))./h(1));
c1 .*=(G(1)./6).*(((xj(3).-x1).^3)./h(2).-(h(2).*(xj(3).-x1))).+(G(2)./6).*(((x1.-xj(2)).^3)./h(2).-(h(2).*(x1.-xj(2)))).+((yj(2).*(xj(3).-x1).+yj(3).*(x1.-xj(2)))./h(2));
c2 .*=(G(2)./6).*(((xj(4).-x2).^3)./h(3).-(h(3).*(xj(4).-x2))).+(G(3)./6).*(((x2.-xj(3)).^3)./h(3).-(h(3).*(x2.-xj(3)))).+((yj(3).*(xj(4).-x2).+yj(4).*(x2.-xj(3)))./h(3));
c3 .*=(G(3)./6).*(((xj(5).-x3).^3)./h(4).-(h(4).*(xj(5).-x3))).+(G(4)./6).*(((x3.-xj(4)).^3)./h(4).-(h(4).*(x3.-xj(4)))).+((yj(4).*(xj(5).-x3).+yj(5).*(x3.-xj(4)))./h(4));
c4 .*=(G(4)./6).*(((xj(6).-x4).^3)./h(5).-(h(5).*(xj(6).-x4))).+(G(5)./6).*(((x4.-xj(5)).^3)./h(5).-(h(5).*(x4.-xj(5)))).+((yj(5).*(xj(6).-x4).+yj(6).*(x4.-xj(5)))./h(5));
c5 .*=(G(5)./6).*(((xj(7).-x5).^3)./h(6).-(h(6).*(xj(7).-x5))).+(G(6)./6).*(((x5.-xj(6)).^3)./h(6).-(h(6).*(x5.-xj(6)))).+((yj(6).*(xj(7).-x5).+yj(7).*(x5.-xj(6)))./h(6));
c6 .*=(G(6)./6).*(((xj(8).-x6).^3)./h(7).-(h(7).*(xj(8).-x6))).+(G(7)./6).*(((x6.-xj(7)).^3)./h(7).-(h(7).*(x6.-xj(7)))).+((yj(7).*(xj(8).-x6).+yj(8).*(x6.-xj(7)))./h(7));
c7 .*=(G(7)./6).*(((xj(9).-x7).^3)./h(8).-(h(8).*(xj(9).-x7))).+(G(8)./6).*(((x7.-xj(8)).^3)./h(8).-(h(8).*(x7.-xj(8)))).+((yj(8).*(xj(9).-x7).+yj(9).*(x7.-xj(8)))./h(8));
c8 .*=(G(8)./6).*(((xj(10).-x8).^3)./h(9).-(h(9).*(xj(10).-x8))).+(G(9)./6).*(((x8.-xj(9)).^3)./h(9).-(h(9).*(x8.-xj(9)))).+((yj(9).*(xj(10).-x8).+yj(10).*(x8.-xj(9)))./h(9));
c9 .*=(G(9)./6).*(((xj(11).-x9).^3)./h(10).-(h(10).*(xj(11).-x9))) .+ ((yj(10).*(xj(11).-x9).+yj(11).*(x9.-xj(10)))./h(10));
y=c;
y1=c1;
y2=c2;
y3=c3;
y4=c4;
y5=c5;
y6=c6;
y7=c7;
y8=c8;
y9=c9;
plot(x,y,x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,x7,y7,x8,y8,x9,y9);
set(gca,'FontSize',20);
xlabel("X",'FontSize',20);
ylabel("g(x)",'FontSize',20);
title("g(x) with 11 datas in hw1CD",'FontSize',20);
grid on;
留言列表
STM32 (1)
