www.isutp.fora.pl

misia

interpolacja zwykla
clear;
clf;
funcprot(0);
x=[-0.57;-0.21;0.11;0.34;0.57];
disp(x,"x = ");
y=[-0.11;0.24;0.34;0.15;-0.36];
disp(y,"y = ")
wm=length(x);
disp("wm = "+string(wm));
wn=wm-1
disp("wn = "+string(wn));
X=ones(wm,wm);
for i=2:5
X(:,i)=x^(i-1);
end
disp(X,"X = ");
c=X\y;
disp(c,"c = ");
function w=wnpolynomial(v);
t=ones(wm,1);
for i=2:wm
t(i)=v^(i-1);
end;
w=sum(c.*t);
endfunction;
d=linspace(x(1)-0.2,x(wm)+0.2,100);
plot(d,wnpolynomial);
plot(x,y,"+k");

interpolacja dla wezlow rownoodleglych
clear;
clf;
y_i=[-1.71;-0.94;-0.11;0.24;0.51;-0.17];
disp(y_i,"y_i = ");
m=length(y_i);
disp(m,"m = ");
n=m-1;
disp(n,"n = ");
h=0.5;
disp(h,"h = ");
x_i=zeros(m,1);
x_i(1)=-1;
for i=2:m
x_i(i)=x_i(i-1)+h;
end;
disp(x_i,"x_i = ");
deltaik=zeros(m,m);
deltaik(:,1)=y_i;
for k=2:m
for i=1 Sad

m-k+1)
deltaik(i,k)=deltaik(i+1,k-1)-deltaik(i,k-1)
end
end;
disp(deltaik,"delta = ");
function om=omega(i, x);
om=x-x_i(1);
for j=2 Sad

i+1);
om=om*(x-x_i(j));
end;
endfunction;
function w=wnpolynomial(x);
w=deltaik(1,1);
for i=0 Sad

n-1)
w=w+(deltaik(1,i+2)*omega(i,x))/(factorial(i+1)*(h^(i+1)))
end
endfunction;
d=linspace(x_i(1)-0.2,x_i(m)+0.2,100);
plot(d,wnpolynomial);
plot(x_i,y_i,"+k");

newton dla dowolnych wezlow
clear;
clf;
x_i=[-1;-0.87;-0.63;-0.32;0.11;0.34;0.72;0.86;1];
f_i=[-1.04;-0.98;-1.1;-0.95;-1.01;-0.9;-0.97;-1;-1.01];
disp(f_i,"f_i = ");
disp(x_i,"x_i = ");
m=length(f_i);
disp(m,"m = ");
n=m-1;
disp(n,"n = ");

fik=zeros(m,m);
fik(:,1)=f_i;
for k=2:m
for i=1 Sad

m-k+1)
fik(i,k)=(fik(i+1,k-1)-fik(i,k-1))/(x_i(i+k-1)-x_i(i));
end
end;
disp(fik,"fik = ");
function om=omega(i, x);
om=x-x_i(1);
for j=2 Sad

i+1);
om=om*(x-x_i(j));
end;
endfunction;
function w=wnpolynomial(x);
w=fik(1,1);
for i=0 Sad

n-1)
w=w+(fik(1,i+2)*omega(i,x))
end
endfunction;
d=linspace(x_i(1)-0.2,x_i(m)+0.2,100);
plot(d,wnpolynomial);
plot(x_i,f_i,"+k");

lagrang dla dowolnych wezlow
clear;
clf;
n=13;
m=n+1;
disp("m = "+string(m));
disp("n = "+string(n));
h=2/n
xi=zeros(m,1);
for i=0:n
xi(i+1)=-1+(i*h)
end
disp(xi,"xi= ")
function f=f(x);
f=1/(1+25*(x^2));
endfunction
fi=zeros(m,1);
for i=1:m
fi(i)=f(xi(i));
end
disp(fi,"fi=");
disp(xi,"xi = ");
fxik=zeros(m,m);//liczymy deltę
fxik(:,1)=fi;
for k=2:m
for i=1 Sad

m-k+1)
fxik(i,k)=(fxik(i+1,k-1)-fxik(i,k-1))/(xi(i+k-1)-xi(i));
end
end
disp(fxik,"fxik = ");
function om=omega(i, x);//obliczamy w
om=x-xi(1);
for j=2 Sad

i+1)
om=om*(x-xi(j));
end
endfunction
function w=wnpolynomial(x);//liczymy wn
w=fxik(1,1);
for i=0 Sad

n-1)
w=w+(fxik(1,i+2)*omega(i,x));
end
endfunction
d=linspace(xi(1),xi(m),100);//rysujemy wykres
plot(d,wnpolynomial);
plot(xi,fi(:,1),"+k");
clear;

czebyszew
clear;
clf;
n=13;
m=n+1;

disp("m = "+string(m));
disp("n = "+string(n));
h=2/n
xi=zeros(m,1);
for i=1:m
xi(i)=-cos(((2*(i-1)+1)/2 Cool

*%pi)
end
disp(xi,"xi= ")
function f=f(x);
f=1/(1+25*(x^2));
endfunction

fi=zeros(m,1);
for i=1:m
fi(i)=f(xi(i));
end
disp(fi,"fi=");
disp(xi,"xi = ");
fxik=zeros(m,m);//liczymy deltę
fxik(:,1)=fi;
for k=2:m
for i=1 Sad

m-k+1)
fxik(i,k)=(fxik(i+1,k-1)-fxik(i,k-1))/(xi(i+k-1)-xi(i));
end
end
disp(fxik,"fxik = ");
function om=omega(i, x);//obliczamy w
om=x-xi(1);
for j=2 Sad

i+1)
om=om*(x-xi(j));
end
endfunction
function w=wnpolynomial(x);//liczymy wn
w=fxik(1,1);
for i=0 Sad

n-1)
w=w+(fxik(1,i+2)*omega(i,x));
end
endfunction
d=linspace(xi(1),xi(m),100);//rysujemy wykres
plot(d,wnpolynomial);
plot(xi,fi,"+r");

clear;

phi1
clear;
clf;
funcprot(0);
a=-1;
b=1;
n=5;
h=(b-a)/n;
function y=phi1(i,x);
xi=zeros(3);
for k=1:3
xi(k)=a+(i-2+k)*h;
end
disp(xi);
if x<xi(1) then
y=0;
else
if x<xi(2) then
y=(x-xi(1))/h;
else
if x<xi(3) then
y=(xi(3)-x)/h;
else
y=0;
end;
end;
end;
endfunction;
d=linspace(a,b,100);
plot(d,phi1(1,x));
clear;

phi3
clear;
clf;
funcprot(0);
a=-1;
b=1;
n=13;
h=(b-a)/n;
function y=phi3(i, x);
for k=1:5
xi(k)=a+(i-3+k)*h;
end
if x<xi(1) then y=0
else
if x<xi(2) then
y=((x-xi(1))^3)/(h^3);
else
if x<xi(3) then
y=1+3*((x-xi(2))/h)+3*(((x-xi(2))^2)/h^2)-3*(((x-xi(2))^3)/h^3);
else
if x<xi(4) then
y=1+(3*(xi(4)-x)/h)+3*(((xi(4)-x)^2)/h^2)-3*(((xi(4)-x)^3)/h^3);;
else
if x<xi(5) then
y=(((xi(5)-x)^3)/h^3);
else
y=0
end
end
end
end
end

endfunction;
dx=linspace(a,b,100);
for j=1:length(dx)
dy(j)=phi3(1,dx(j));
end
plot(dx,dy);

scf(1);
for i=-1:n+1;
dx=linspace(a,b,100);
for j=1:length(dx);
dy(j)=phi3(i,dx(j));
end
plot(dx,dy);
end
clear;