%interpolowanie: Lagrange: function L=Lagrange (x,y,p) %obl wart. wiel. Lagrange'a dla pod. arg p
n=length(x)
L=0;
il=1;
for i=1:1:n %dotyczy sumy
for k=1:1:n %dotyczy iloczynu
if k~=i %pętla warunkowa
il=il*(p-x(k))/(x(i)-x(k));
end;
end;
L=L+y(i)*il;
end;
end
%Całkowanie numeryczne: TP,SP,TZ,SZ:
clc;
a=pi; %poczatek przed calkowania
b=2*pi; %koniec przed calkowania
f=inline('sin(x)/x'); %nasza funkcja
T=(b-a)/2*(f(a)+f(b)); %met trapezow
S=(b-a)/6*(f(a)+4*f((a+b)/2)+f(b)); %met simsona
disp(['Wartosc calki wynosi'])
disp(['dla metody trapezow prostej: ', num2str(T)])
disp(['dla metody Simsona prostej: ', num2str(S)])
m=8; %liczba czesci
h=(b-a)/m; %dlugosc przedzialu
sum=0;
for i=1:1:m-1; %petla obliczajaca wartosc w punktach pomiedzy a,b
sum=sum+f(a+i*h);
end;
Tz=h/2*(f(a)+f(b)+2*sum); %wzor na met zlozona simsona
disp(['dla metody trapezow zlozonej z ', num2str(m), ' czesci: ', num2str(Tz)])
sum1=0;
for i=1:1:m-1;
sum1=sum1+f(a+i*h);
end;
sum2=0;
for i=1:1:m;
sum2=sum2+f(a+(i-0.5)*h);
end;
Sz=h/6*(f(a)+f(b)+2*sum1+4*sum2);
disp(['dla metody simsona zlozonej z ', num2str(m), ' czesci: ', num2str(Sz)])
%Rozw.równ z 1niewiadom: Bisekcja pierwiastek(brak Newto):
clc;
format long;
f=inline('cos(x)-2*x');
a=0;
b=1;
k=0;
maxn=100;
eps=10e-10;
ezplot(f,[-10,10])
grid;
if f(a)*f(b)>0;
disp(['Blad - funkcja f nie ma pierwiastów na tym przedziale']);
return %lub break
end;
while k<maxn && abs(b-a)>=eps;
k=k+1;
x=(a+b)/2;
if f(a)*f(x)<=0;
b=x;
else a=x;
end;
disp(['Pierwiastkiem tej funkcji dla k=', Num2Str(k) ' jest liczba ']);
x
disp(['F(x)=']);
f(x)
end;
%Rozw. Równ.różn:1)Euler:
clc;
clear;
rd=inline('exp(3/2*t^2)')
hold on;
grid;
ezplot(rd,[0,1.5]);
plot(0,1,'ro');
f=inline('3*t*y');
a=0;
b=1.5;
y(1)=1
n=10;
h=(b-a)/n;
t=a:h:b;
for i=2:1:length(t);
y(i)=y(i-1)+h*f(t(i-1),y(i-1))
plot(t(i),y(i),'ro');
end;
%2)RK:
clc; clear;
rd=inline('exp(3/2*t^2)')
hold on;
grid;
ezplot(rd,[0,1.5,0,40]);
plot(0,1,'ro');
f=inline('3*t*y');
a=0;
b=1.5;
y(1)=1
n=10;
h=(b-a)/n;
t=a:h:b;
for i=2:1:length(t);
K1=h*f(t(i-1),y(i-1))
K2=h*f(t(i-1)+h/2,y(i-1)+K1/2)
K3=h*f(t(i-1)+h/2,y(i-1)+K2/2)
K4=h*f(t(i-1)+h,y(i-1)+K3)
y(i)=y(i-1)+(K1+2*K2+2*K3+K4)/6
plot(t(i),y(i),'ro');
end;