www.isutp.fora.pl

misia

//rownoodlegle
clear;
clf;
funcprot(0);
a=-1;
b=1;
n=11;
al1=1.634;
bt1=-4.673;
h=(b-a)/n;
xi=zeros(n+1,1);
for i=0:n
xi(i+1)=a+i*h;
end;
disp(xi,"xi = ");
function f=f(x);
f=1/(1+25*(x^2));
endfunction
yi=zeros(n+3,1);
yi(1)=-(h*al1)/3;
for i=0:n
yi(i+2)=f(x(i+1));
end;
yi(n+3)=-(h*bt1)/3;
disp(yi,"yi = ");
M=zeros(n+3,n+3);
M(1,1)=1;
M(1,3)=-1;
for i+2 Sad

n+2)
for j=1 Sad

n+3)
if i=j then M(i,j)=4;
end;
if i+1=j|i=j+1 then M(i,j)=1
end;
end;
end;
M(n+3,n+1)=1;
M(n+3,n+3)=-1;
disp(M,"M= ");

2 sposob

clear;
clf;
funcprot(0);
a=-1;
b=1;
n=13;
h=(b-a)/n;
al2=0.1;
bt2=-0.1;
xi=zeros(n+1,1);
for i=0:n
xi(i+1)=a+i*h;
end;
disp(xi,"xi = ");
function f=f(xi);
f=1/(1+25*(xi^2));
endfunction;
yi=zeros(n+3,1);
yi(1)=-(h*al2)/6;
for i=0:n
yi(i+2)=f(xi(istyle="color: #5c5c5c;">+1));
end
yi(n+3)=-(h*bt2)/6;

disp(yi,"yi = ");

M=zeros(n+3,n+3);
M(1,1)=1;
M(1,2)=-2;
M(1,3)=1;
for i=2 Sad

n+2)
for j=1 Sad

n+3)
if i==j then M(i,j)=4;
end
if i+1==j|i==j+1 then M(i,j)=1;
end;

end;
end;
M(n+3,n+1)=1;
M(n+3,n+2)=-2;
M(n+3,n+3)=1;

disp(M,"M = ");

c=M\yi;
disp(c,"c = ");

function y=phi3(i, x);
for k=1:5
xi(k)=a+(i-3+style="color: #000000;">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*(((style="color: #000000;">xi(4)-x)^2)/h^2)-3*(((xi(4)-x)^3)/h^3);;
else
if x<xi(5) then
y=(((xi(5)-xstyle="color: #4a55db;">)^3)/h^3);
else
y=0
end
end
end
end
end

endfunction;
function s=s(x);
s=0;
for i=-1 Sad

n+1);
s=s+c(i+2)*style="color: #000000; text-decoration: underline;">phi3(i,x);
end;
endfunction;

d=linspace(a,b,100);
plot(d,f,"k");
plot(d,s,"r");
plot(xi,yi(2 Sad

n+2style="color: #4a55db;">)),"+o");

//sklejanie dowolne wezly
clear;
clf;
funcprot(0);
a=-1;
b=1;
n=13;
m=n+1;
xi=zeros(m,1);
xi(1)=a;
xi(m)=b;

function tr=tr(t);
tr=((b-a)*t+(a+b))/2;
endfunction;
for i=0 Sad

n-2)
xi(i+2)=tr(-cos(((2*i+1)*%pi)/(2*(n-1))));
end;
function f=f(xi);
f=1/(1+25*(xi^2));
endfunction;
yi=zeros(m,1);
for i=1:m
yi(i)=f(xi(i));
end;

disp(xi,"xi = ");
disp(yi,"yi = ");

// twierdzenie 2.6,
hi=zeros(n);
for i=1:n
hi(i)=xi(i+1)-xi(i);
end
disp(hi,"hi = ");
di=zeros(n);
for i=1 Sad

n-1)
di(i)=(6/(hi(i)+hi(i+1)))*((yi(i+2)-yi(i+1))/hi(i+1)-((yi(i+1)-yi(i))/hi(i)));
end;
di(n)=(6/(hi(1)+hi(n)))*((yi(2)-yi(1))/hi(1)-((yi(n+1)-yi(n))/hi(n)));
disp(di,"di = ");

lai=zeros(n);
for i=1 Sad

n-1)
lai(i)=(hi(i+1))/(hi(i)+hi(i+1));
end;
lai(n)=(hi(1))/(hi(1)+hi(n));
disp(lai,"lai = ");

mi=zeros(n);
for i=i:n
mi(i)=1-lai(i);
end;
disp(mi,"mi = ");

mg=zeros(n+1,n+1);
mg(1,1)=1;
mg(1,n+1)=-1;
for i=1 Sad

n-1)
mg(i+1,i)=mi(i);
mg(i+1,i+1)=2;
mg(i+1,i+2)=lai(i);
end
mg(n+1,1)=2*lai(n);
mg(n+1,2)=lai(n);
mg(n+1,n)=mi(n);
mg(n+1,n+1)=2*mi(n);
disp(mg,"mg = ");

ww=zeros(n+1,1); //wktor wyrazow wolnych
ww(1)=0;
for i=1 Sad

n-1)
ww(i+1)=di(i);
end
ww(n+1)=di(n);
disp(ww,"ww = ");

M=mg\ww;
disp(M,"M=");

// Aproksymowac funkcje f(x)=x^sinx na przedziale <0,2> wykorzystujac do tego wielomiany legendra na przedziale <-1,1>

// 1) zadeklarowac f
// 2) wprowadzic a i b
// 3) wybrac i wprowadzic m, N
// m>=0 N>=2
// 4)zadeklarowac funkcje tr i tr^-1
// 5)zadeklarowac g(t)=f(tr^-1)
// G(t)=sum od j=1 do m c_jP_j(t)
//legendre(i,0,t)
//P_i(t)
//6) założyc wektor c odpowiedniego wymiaru
//7) założyc wektor tv=linspece(-1,1,N)
//C_i=(f,P_i)/(P_i,p_i)
// Cool

dla j=0,...,m obliczam (g,P_j) i (P_j,P_j) i C_i=(g,P_j)/(P_j,p_j)
//założyc wektor gP z tym wymiarem co wektor z 7)
//gP(i)=g(tr(i))*(legendre(j,0,tr(i)))^2

clear;
funcprot(0);
m=2;
disp("m = "+string(m));
N=25;
disp("N = "+string(N));

function tr=tr(x);
tr=(2*x-(a+b))/(b-a)
endfunction;

function tri=tr_inv(t);
tri=((b-a)*t+(a+b))/2;
endfunction;

function g=g(t);
g=f(tri(t))
endfunction;

f12i(i)=f1(t_i(i))*f2(t_i(i));
c=zeros(m+1);
l=nearfloat("pred",1)
t_i=linspace(-1,1,N);
lci=zeros(N);
mci=zeros(N);
for k=0:m
for i=1:N
lci(i)=g(t_i(i))*legendre(k,0,t_i(i))
mci(i)=g(t_i(i))*legendre(k,0,t_i(i))

end
end
c(k+1)=inttrap(t_i,lci)
function G(t)
if t<= then t=nearfloat("succ",-1)
if t>=1 then t=nearfloat("pred",1)
end
end
g=sum(cp)
endfunction

//całka od -1 do 1 h(t)dt
//funkcja tworzę wektor t=linspace (-1,1,N)
//utwórz drugi wektor ht z warunkiem ht (ht(i)=h(t(i)))
//funkcja całki intrap(t,ht)

//aproksymacja dowolnymi węzłami
clear;
funcprot(0);
a=0;
disp("a="+string(a));
b=2;
disp("b="+string(b));
m=1;
disp("m="+string(m));
N=20;
function f=f(x);
f=x^(sin(x));
endfunction;
function tri=tri(t);
tri=((b-a)*t+(a+b))/2;
endfunction;
function tr=tr(x);
tr=(2*x-(a+b))/(b-a);
endfunction;
function g=g(t);
g=f(tri(t));
endfunction;
L=nearfloat("pred",1);
tc=linspace(-L,L,N);
c=zeros(m+1);
lc=zeros(N);
mc=zeros(N);
for k=0:m
for i=1:N
lc(i)=g(tc(i))*legendre(k,0,tc(i));
mc(i)=legendre(k,0,tc(i))^2;
end;
c(k+1)=inttrap(tc,lc)/inttrap(tc,mc);
end;
function G=G(t);
//if t<=-1 then
//t=nearfloat("succ",-1);
//if t>=1 then
//t=nearfloat("pred",-1);
G=0
for k=0:m
G=G+(c(k+1)*legendre(k,0,t));
end
endfunction;
function F=F(x);
F=G(tr(x));
endfunction;
d=linspace(a,b,100);
plot(d,f);
plot(d,F,"r");
disp(c,"c=");

// dyskretna sinusowa transformacja

clear;
clf;
funcprot(0);
a=0;
b=2;
c=1;
N=15;

disp(a,"a = ");
disp(b,"b = ");
disp(c,"c = ");
disp(N,"N = ");

function f=f(x)
if x<0 then f=0; else
if x<(1/2) then f=1; else
if x<=1 then f=3; else
if x<(3/2) then f=-2*x+5; else
if x==(3/2) then f=0; else
if x<2 then f=2; else
if x==2 then f=1; else
f=0;
end
end
end
end
end
end
end
endfunction;

function tr=tr(x);
tr=(%pi/(b-a))*(x-a);
endfunction

function tri=tri(t);
tri=(((b-a)/%pi)*t)+a;
endfunction

function g=g(t);
g=f(tri(t))-c;
endfunction

xk=zeros(N);
gk=zeros(N);
Gn=zeros(N);
bn=zeros(N);

for i=0 Sad

N-1)
xk(i+1)=(%pi*i)/N;
gk(i+1)=g(xk(i+1));
end;

disp(xk,"xk = ");
disp(gk,"gk = ");

for j=0 Sad

N-1)
gks=zeros(N-1);
for i=1 Sad

N-1)
gks(i)=gk(i+1)*sin(xk(i+1)*j);
end;
Gn(j+1)=sum(gks);
bn(j+1)=(2/N)*Gn(j+1);
end;

disp(Gn,"Gn = ");
disp(bn,"bn = ");

function G=G(t);
bns=zeros(N-1);
for j=1 Sad

N-1)
bns(j)=bn(j+1)*sin(j*t);
end
G=sum(bns);
endfunction;

function F=F(x);
F=G(tr(x))+c;
endfunction;

function E=E(x);
E=abs(f(x)-F(x));
endfunction;

h=b-a;
disp(h,"h = ");
dl=linspace(a-2*h,b+2*h,5000);
d=linspace(a,b,1000);
d1=linspace(a,0.5-%eps,250);
d2=linspace(0.5,b-%eps,750);

scf(0);
//plot(d,f);
plot(d1,f,);
plot(d2,f);
plot(dl,F,"r");

plot({0,0.5,1,1.5,2},{1,3,3,0,1},".");
plot({0.5,1.5,2},{1,2,2},"o");

scf(1);
plot(d,F,"r");
//plot(d,f,"b");
plot(d1,f);
plot(d2,f);
plot({0,0.5,1,1.5,2},{1,3,3,0,1},".");
plot({0.5,1.5,2},{1,2,2},"o");

scf(2);
//plot(d,E,"g");

plot(d1,E,"r");
plot(d2,E,"r");
plot({0.5,1.5,2},{abs(1-F(0.5)),abs(2-F(1.5)),abs(2-F(2))},"o")
plot({0.5,1.5,2},{E(0.5),E(1.5),E(2)},".");