[ Pobierz całość w formacie PDF ]
% Wielomian interpolacyjny metoda Lagrangeaclearclchold ongrida=0b=6%axis([a,b,a,b])f=inline('2*sin(pi/6*x)')x=[0,1,3,6]y=f(x)plot(x,y,'ro')x2=a:0.1:b;y2=f(x2);plot(x2,y2,'g-') %rysuje f-cjeyw=Lagrange(x2)%n=length(x) %dｳugo懈 wektora, dｳuｿszy wymiar macierzy%w=polyfit(x,y,n-1) %wyznacza w黝ｳy interpolacji%yw=polyval(w,x2) %oblicza warto懈 wielomianu W dla podanego argumentuplot(x2,yw,'r-') %rysuje wielomian interpolacyjnyp=input('Podaj p= ') %wprowadzanie zmiennych z klawiatury%yp=polyval(w,p)yp=Lagrange(x,y,p)disp('Warto懈 wielomianu dla podanego argumentu wynosi ')disp(yp)% Wielomian interpolacyjny:clccleargrid;a=-1; b=1;n=3; %3,5,11,15 %liczba w黝ｳf=inline('abs(x)');%WﾊZ｣Y RﾓWNOODLEG｣E:x=a:(b-a)/(n-1):b %w黝ｳyy=f(x)wr = polyfit(x,y,n-1) %wyzn. wspcz. wiel. interp.xrys = a:0.01:b;yrysR = polyval(wr,xrys); %obl. wart. wiel.%WﾊZ｣Y CZEBYSZEWA:for k=0:1:n-1c(k+1)=(b-a)/2*cos((2*k+1)/(2*n)*pi)+(a+b)/2endyc=f(c)wc = polyfit(c,yc,n-1) %wyzn. wspcz. wiel. interp.yrysC = polyval(wc,xrys); %obl. wart. wiel.plot(xrys,f(xrys),'r-',x,y,'go',xrys,yrysR,'g-',c,yc,'bo',xrys,yrysC,'b-')p=0.95;disp(['wartosci dla argumentu =', num2str(p)])disp(['f(x)=', num2str(f(p))])disp(['wr(x)=',num2str(polyval(wr,p))])disp(['wc(x)=',num2str(polyval(wc,p))])% Caｳki (metoda trapez I Simsona)clc;a=pi; %poczatek przed calkowaniab=2*pi; %koniec przed calkowaniaf=inline('sin(x)./x'); %nasza funkcjaT=(b-a)/2*(f(a)+f(b)); %met trapezowS=(b-a)/6*(f(a)+4*f((a+b)/2)+f(b)); %met simsonadisp(['Wartosc calki wynosi'])disp(['dla metody trapezow prostej: ', num2str(T)])disp(['dla metody Simsona prostej: ', num2str(S)])m=3; %liczba czescih=(b-a)/m; %dlugosc przedzialusum=0;for i=1:1:m-1; %petla obliczajaca wartosc w punktach pomiedzy a,bsum=sum+f(a+i*h);end;Tz=h/2*(f(a)+f(b)+2*sum); %wzor na met zlozona simsonadisp(['dla metody trapezow zlozonej z ', num2str(m), ' czesci: ', num2str(Tz)])sum1=0;for i=1:1:m-1; %petla obliczajaca wartosc w punktach pomiedzy a,bsum1=sum1+f(a+i*h);end;sum2=0;for i=1:1:m; %petla obliczajaca wartosc w punktach pomiedzy a,bsum1=sum1+f(a+i*h);end;% Trapez zｳoｿony, Simson zｳoｿony:clc;a=pi; %poczatek przed calkowaniab=2*pi; %koniec przed calkowaniaf=inline('1/sqrt(x.^3-1)'); %nasza funkcjaT=(b-a)/2*(f(a)+f(b)); %met trapezowS=(b-a)/6*(f(a)+4*f((a+b)/2)+f(b)); %met simsonadisp(['Wartosc calki wynosi'])disp(['dla metody trapezow prostej: ', num2str(T)])disp(['dla metody Simsona prostej: ', num2str(S)])m=32; %liczba czescih=(b-a)/m; %dlugosc przedzialusum=0;for i=1:1:m-1; %petla obliczajaca wartosc w punktach pomiedzy a,bsum=sum+f(a+i*h);end;Tz=h/2*(f(a)+f(b)+2*sum); %wzor na met zlozona simsonadisp(['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)])%Caｳkowanie numeryczne: TP,SP,TZ,SZ:clc;a=pi; %poczatek przed calkowaniab=2*pi; %koniec przed calkowaniaf=inline('sin(x)/x'); %nasza funkcjaT=(b-a)/2*(f(a)+f(b)); %met trapezowS=(b-a)/6*(f(a)+4*f((a+b)/2)+f(b)); %met simsonadisp(['Wartosc calki wynosi'])disp(['dla metody trapezow prostej: ', num2str(T)])disp(['dla metody Simsona prostej: ', num2str(S)])m=8; %liczba czescih=(b-a)/m; %dlugosc przedzialusum=0;for i=1:1:m-1; %petla obliczajaca wartosc w punktach pomiedzy a,bsum=sum+f(a+i*h);end;Tz=h/2*(f(a)+f(b)+2*sum); %wzor na met zlozona simsonadisp(['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n 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 na tym przedziale']);return %lub breakend;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 ']);xdisp(['F(x)=']);f(x)end;%Rozw. Rn.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)=1n=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)=1n=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)/6plot(t(i),y(i),'ro');end; [ Pobierz całość w formacie PDF ]