Metodat e Analizës Numerike në Matlab Versioni 1
MENTOR SHEVROJA
Tirane 2013
Permbajtja Hyrje .............................................................................................................................................................. 4 1.
Njohuri paraprake ................................................................................................................................. 5 1.1
Variablat ........................................................................................................................................ 5
1.2
Vektoret ........................................................................................................................................ 5
1.3
Matricat......................................................................................................................................... 6
1.4
Struktura e kontrollit..................................................................................................................... 7
1.5
Strukturat ciklike ........................................................................................................................... 8
1.6
Funksionet..................................................................................................................................... 9
1.7
Skripte ......................................................................................................................................... 10
1.8
Komanda ndihmese .................................................................................................................... 12
2.
Gabimi absolut dhe relativ .................................................................................................................. 13
3.
Zgjidhja numerike e ekuacioneve algjebrike dhe transhengjente ...................................................... 14
4.
5.
6.
3.1
Metoda e pergjysmimit............................................................................................................... 14
3.2
Metoda e Njutonit ...................................................................................................................... 16
3.3
Metoda e kordes (Lagranzhit) ..................................................................................................... 17
3.4
Metoda e pikes fikse ................................................................................................................... 17
Sisteme te ekuacioneve lineare .......................................................................................................... 18 4.1
Zgjidhja e sistemeve me matrice trekendeshe te poshtme ........................................................ 18
4.2
Metoda e Gausit.......................................................................................................................... 18
4.3
Faktorizimi LU ............................................................................................................................. 20
4.4
Metoda e Jakobit ........................................................................................................................ 20
Interpolimi i funksioneve .................................................................................................................... 22 5.1
Polinomi interpolues i Lagranzhit ............................................................................................... 22
5.2
Polinomi interpolues i Njutonit ................................................................................................... 23
Perafrimi i funksioneve ....................................................................................................................... 25 6.1
7.
8.
Metoda e katroreve me te vegjel ............................................................................................... 25
Integrimi i perafert i funksioneve ....................................................................................................... 28 7.1
Metoda e perbere e trapezit....................................................................................................... 28
7.2
Metoda e perbere e Simpsonit ................................................................................................... 29
Zgjidhja e perafert e ekuacioneve diferenciale te zakonshme ........................................................... 31 8.1
Metoda e Eulerit ......................................................................................................................... 31 2
8.2
Metoda Runge-Kutta................................................................................................................... 32
3
Hyrje Matlab-i eshte nje program te nje niveli te lart ne fushen e njehsimit numerik, me shume funksione te gatshme qe ofrojne lehtesi ne zgjidhjen e metodave numerike. Nepermjet Matlabit mund te zhvillojme algoritme, te analizojme te dhena, te krijojme module dhe aplikime te ndryshme. Ne kete material do fokusohemi ne pjesen e zhvillimit te algoritmeve per Metodave te Analizës Numerike.
4
1. Njohuri paraprake 1.1
Variablat
Ne ambientin e punes se Matlab-it te gjitha variablat jane matrica, dhe tipi i te dhenave te zbatuara (default) eshte me precision te dyfishte (double). Çdo gje qe ne deklarojme ne Matlab eshte matrice, me kete lidhet edhe emri i programit Matlab (Matrix Laboratory). Matricat me permasa 1x1 do ti quajme variabla skalar ose thjeshte variabel. Deklarimi i variablave ne Matlab eshte shume i thjeshte, mjafton te shkruajme >> x=1.75 x= 1.7500 Ne kete menyre kemi deklaruar nje variabel dhe i kemi dhene vleren 1.75. Po te vini re, pasi kemi shkruar variabin dhe i kemi dhene vleren 1.75, Matlab-i ka afishuar automatikisht emrin e variablit dhe vleren e tij. Nese deshirojme te afishojme vleren e nje variabli atehere mjafton te shkruajme emrin e tij ne dritaren e komandave, Matlab-i automatikisht do te afishoj emrin e variablit dhe vleren e tij. Rregulli per te emertuar nje variabel eshte te filloj me nje shkronje, vazhdimi mund te jete shkronje ose numer ose underscore (_) dhe me gjatesi jo me te madhe se vlera e konstantes namelengthmax te matlabit. Matlab-i eshte i ndjeshem ndaj shkronjave te medha dhe te vogla (case sensitive), kjo do te thote se x dhe X jane dy variabla te ndryshem. Nese ne perfundim te shprehjes qe shkruani ne dritaren e komandave (Command Window) nuk shenoni pikepresje (;), Matlabi afishon vleren e shprehjes ne dritaren e komandave. Matlabi ka nje variabel special ans qe eshte shkurtimi i answer qe nenkupton pergjigjen e shprehjes se fundit qe ju keni shkruar ne dritaren e komandave te Matlab-it.
1.2
Vektoret
Nje vektor eshte nje varg elementesh nje dimensional. Vektoret mund ti ndajme ne dy tipe: vektor te tipit rresht dhe vektor te tipit shtylle. Vektor rresht shkruhet horizontalisht; vektor shtyll shkruhet vertikalisht. Ne Matlab, vektoret deklarohen duke shkruar elementet ne kllapa katrore. Per vektoret rresht elementet e vektorit shkruhen te ndara me hapesire ose presje, >> x=[5 3 -2 4 -1 0 2] x= 5
3 -2
4 -1
0
2
5
Ndersa vektori shtylle shkruhet ne kllapa katrore por elementet ndahen me pikepresje ndermjet tyre >> x=[5; 3; -2; 4; -1; 0; 2] x= 5 3 -2 4 -1 0 2
Po te shikoni edhe afishimi i vektoreve behet ne varesi se si eshte deklaruar vektori rresht apo shtylle. Nese nuk duam te na shfaqet vektori mjafton qe ne fund te shenojme ‘;’ >> x=[5 3 -2 4 -1 0 2]; Kur duam te aksesojme nje element te caktuar te vektorit, atehere shkruajm emrin e vektorit te cilit duam ti aksesojme nje komponente te tij, dhe me pas indeksin e elementit qe duam te aksesojme ne kllapa rrethore (indeksi). Psh ne vektorin e mesiperm x duam te perdorim elementin e trete, shkruajme >> x(3) ans = -2 Sic e shikoni vlera e elementit te trete eshte -2.
1.3
Matricat
Matricat jane koleksione dy dimensionale elementesh te te njejtit tip te dhenash. Deklarimi i matricave behet ne ngjashmeri me vektoret, elementet vendosen ne kllapat katrore duke ndare elementet me hapesire ose presje dhe rreshtat ndahen nga pikepresje (;) >> C=[-2 5 1;1 0 4] C=
6
-2
5
1
1
0
4
Ne kete menyre kemi deklaruar nje matrice C me permasa 2x3. Per te aksesuar elementin ne rreshtin e pare dhe shtyllen e dyte shkruajme >> C(1,2) ans = 5 Gjithashtu mund te aksesojme elementet e nje rreshti ose shtylle duke shkruar ‘:’ ne pozicionin e rreshtit ose te shtylles si psh >> C(:, 3) % afisho shtyllen e trete te matrices C >> C(2, :) % afisho rreshtin e dyte te matrices C Nepermjet ketyre shkurtimeve mund te ndertojme vektore si psh >> x=C(2, :) x= 1
0
4
Me rreshtin e dyte ndertuam nje vektor rrjesht.
1.4
Struktura e kontrollit
Komanda themelore permes se ciles realizohet degezimi ne Matlab eshte komanda if. Gjate shkruarjes se komandes mund te paraqiten disa raste te degezimeve: te zakonshme, te perbera, te shumefishta, te nderthurura. Sintaksa pergjithshme e kesaj komande eshte: if shprehje komanda elseif shphreje komanda else komanda 7
end Ne vend te shprehje mund te vendoset nje shprehje aritmetiko-logjike ndersa ne vend te komanda mund te vendoset nje ose disa komanda te Matlab-it. Ne sintaksen e pergjithshme te komandes if mund te mungoj njeri ose te dy opsionet elseif dhe else. Kjo komande zbatohet kryesisht kur hasen probleme qe permbajne kerkesen ne qofte se ... atehere... per ndryshe ... . Nje shembull konkret per kete komande me te tre degezime eshte x=input('vendos nje numer x = '); y=input('vendos nje numer tjeter y = '); if x > y disp('x > y') elseif x < y disp('x < y') else disp('x = y') end
Ne fillim perdoruesi vendos dy numra dhe programi pasi gjen raportin e ketyre dy numrave midis tyre e afishon per perdoruesin. Komandat e mesiperme jane shkruar ne nje skedare. Me poshte do te shohim se si krijohen keto skedare me komanda te Matlab-it qe do ti quajme skedare me komanda (script file).
1.5
Strukturat ciklike
Cikli for Cikli for perdoret per te ekzekutuar nje numer te caktuar heresh nje bllok me komanda. Te kryhen 5 iteracione per shprehjen
x=3; for i=1:5 x=x/2+1/x; disp(x); end
Ne kete cikel for-i do te perseritet 5 here. Heren e pare i do kete vleren 1, heren e dyte 2 e keshtu me radhe derisa te behet 5 heren e fundit.
Cikli for perdoret kryesith per te bredhur neper elementet e vektorit ose matrices, ku numri hereve qe ekzekutohet kodi dihet paraprakisht.
8
Ne matricen
zbrisni nga rreshti i dyte rreshtin e pare
A=[1 5 3 5; 6 4 -2 5;6 0 -5 2]; [m,n]=size(A); % m = numri i rreshtave, n = numri i shtyllave. for i=1:n A(2,i)=A(2,i)-A(1,i); end disp(A) % afisho matricen
Cikli while Ky cikel perdoret per te perseritur ekzekutimin e nje kodi ku numri i perseritjeve percaktohet nga rezultatet e llogaritjeve ne trupin e vet ciklit. x=input('Vendos nje numer pozitiv x '); m=0.8; while x>=1 x=x*m; disp(x); end disp('fund');
Cikli i fundit while varet nga numri qe perdoruesi vendos. Sa me i madh te jete numri qe vendos perdoruesi aq me shume iteracione do te kryhen derisa x te behet me i vogel se 1. Me ciklet for dhe while ngusht lidhet edhe komanda break. Kjo komande ben daljen e menjehershme nga cikli.
1.6
Funksionet
Funksioni eshte nje pjese kodi me komanda qe mund ta perdorim disa here. Forma e pergjithshme e funksonit ne Matlab deklarohet ne kete form: function [parametrat dales] = emri_funksionit (parametrat hyres) % eshte mire qe ne krye te shenoni komente se cfare bene ky kod Komandat qe llogarisin parametrat dales. Ne llogartije hyne edhe parametrat hyres nese funksioni ka end % funksioni perfundon me end
Parametrat dales perfaqesojne vlerat qe duam te na kthej funksioni. Parametrat hyres jane vlera qe jane te disponueshme per funksionin qe te llogarisi parametrat dales. Nese ka me shume se nje parameter dales, ato shkruhen ne kllapa katrore [para1 para2]. 9
Shembull. Te shkruajme nje funksion qe te na llogarisi vellimin dhe siperfaqen e nje cilindri me diameter d dhe lartesi h. Informacioni qe funksioni do te perdori jane diametri d dhe lartesia h. Keto parametra do jene te nevojshme per te njehsuar vellimin dhe siperfaqen e cilindrit, keshtu qe do jene parametra hyres per funksionin. Si parametra dales ne kerkojme qe te na njehsoj vellimin v dhe siperfaqen s te cilindrit. function [v, s]=cilinder(d,h) % llogaritja e vellimit v dhe e siperfaqes s te cilindrit % cilindri perfaqesohet nga diametri d dhe lartesi h % rezja sa gjysma e diametrit r=d/2; % velimi eshte siperfaqja e bazes x lartesi v=pi*r^2*h; % siperaqja e cilindrit eshte dyfishi i siperfaqes % se bazes here siperfaqja anesore s=2*pi*r^2+2*pi*r*h; end
Eshte e rendesishme qe kete funksion qe kemi shkruar ta ruajme me te njejten emer si emrin e funksionit. Pasi ta ruajme funksionin jemi gati te ekzekutojme fuksionin dhe te marrim rezultatet ne dritaren e komandave. >> [vellimi, siperfaqja]=cilinder(2,3) vellimi = 9.4248 siperfaqja = 25.1327 >> cilinder(2,3) ans = 9.4248 Sic mund te shihni nga shembulli, kemi disa raste se si mund te marrim rezultatet. I pari eshte te perftojme te gjitha parametra dales si vellimi dhe siperfaqja. Nese nuk percaktojme parametrat dales atehere Matlabi na kthen vleren e parametrit te pare dales si pergjigje te funksionit.
1.7
Skripte
Skripte jane skedare me komanda te Matlab-it. Si komanda te skriptit mund te sherbej jo vetem cdo komande e Matlab-it por edhe funksione te ndertuara prej nesh. Skriptet jane te pershtatshme kur ne duam te ekzekutojme disa komanda njera pas tjetres disa here. 10
Per te krijuar nje skript ne Matlab shkojme tek File > New > Script ose mund te perdorim komanden e shpejte ctrl + n nga tastiera. Menjehere na shfaqet editori ne te cilen shkruajme komandat njera pas tjeteres si ne figure
Pasi te kemi perfunduar shkrimin e komandave ruajme skedarin me nje emer qe duhet ta mbajme mend pasi vetem nepermjet emrit do mund te ekzekutojme komandat qe ndodhet brenda tij. Ne figuren e mesiperme skedari eshte ruajtur me emrin skript1.m. Per ta ekzekutuar kete skedar me komanda, shkruajme emrin e skedarit ne dritaren e komandave dhe perftojme rezultatet. Ne rastin e mesiperm shkruajme skript1, nuk eshte e nevojshme edhe prapashtesa .m.
11
1.8
Komanda ndihmese
Komanda clear Ne situata te ndryshme kur kemi krijuar shume variabla ne hapsiren e punes ne Matlab, eshte e nevojshme te pastrojme te gjitha variablat nga hapsira e punes. Variablat qe ne kemi krijuar afishohen ne dritaren workspace
Po te shkruani ne dritaren e komandave komanden clear do te shikoni qe dritarja e hapesires se punes do jete pastruar nga variablat. Ky pastrim i hapesires se punes eshte i nevojshem ne shume raste per mos pasur ngateresa te variablave kur ne perfundojme pune me nje kod dhe ekzekutojme nje tjeter. Komanda clc Komanda clc ka nje funksion tjeter per te na lehtesuar punen ne Matlab. Kjo komande kryen pastrimin e dritares se komandave. Pasi ne kemi shkruar komanden clc ne dritaren e komandave, Matlab-i pastron dritaren e komandave dhe e kthen ne gjendjen fillestare. Komanda whos Nese shkruani whos ne dritaren e komandave atehere do tju afishohen te gjitha variabla qe ju keni krijuar ne ambientin e punes. Nese komandes whos i percaktoni nje variabel, atehere do tju afishoj informacion rreth variablit. Konsante Ne Matlab ka disa konstante. Keto konstante mund te ndryshohen nese ju perdorni te njejtin emertim te variablit si edhe konstantet. Prandaj duhet te tregoni kujdes qe mos te krijoni variabla me keto emra si pi, eps, i, j. Nese ju krijoni variabel me keto emra atehere kuptimi i vertet i ketyre variablave do te ndryshoj me vleren qe ju keni dhene.
12
2. Gabimi absolut dhe relativ Te ndertohet nje funksion qe te njehsoj gabimin absolut dhe relativ kur jepet numri i sakte dhe vlera e perafert e tij. function [abs_error, rel_error] = errors(x, x0) % funksioni per njehsimin e gabimit absolut dhe relativ % hyrjet: % x-vlera e sakte % x0 - vlera e perafert % daljet: % abs_error - gabimi absolut % rel_error - gabimi relativ abs_error = abs(x-x0); rel_error = abs_error / abs(x); end
Funksion i Matlab-it abs(x-x0)
Njehsimi i vleres absolute te nje numri. Ne kete rast njesohet vlera absolute e diferences se dy numrave x dhe x0
Te modifikohet funksioni i mesiperm per rastin kur duam gabimin relativ te shprehur ne perqindje.
13
3. Zgjidhja numerike e ekuacioneve algjebrike dhe transhengjente 3.1
Metoda e pergjysmimit
Te ndertohet algoritmi per metoden e pergjysmimit ne formen e nje skripti. % % % % % % % %
skripti per metoden e pergjysmimit zgjidhja f(x) = 0 duke perdorur metoden e pergjysmimit Hyrjet: f = funksioni i cili do ti gjendet rrenja a = skaji i majt i intervalit ku ndodhet rrenja b = skaji i djatht i intervalit ku ndodhet rrenja tol = toleranca per ndalimin e procesit iterativ
% Daljet: % x = rrenja e ekucaionit % eps = gabimi me te cilin eshte gjetur rrenja f=input('f='); a=input('a='); b=input('b='); tol=input('tol='); f = inline(f); it=0; fa=f(a); fb=f(b); % % kontrollo nese ne skaje ka shenja te kunderta % if sign(fa)==sign(fb), error('Rrenja nuk ndodhet ne intervalin e dhene'); end % % fillo ciklin % while (b-a)>tol it=it+1; c=(b+a)/2; fc=f(c); if sign(fa)~=sign(fc); % rrenja ndodhet ne [a,c] b=c; fb=fc; else % rrenja ndodhet ne [c,b] a=c; fa=fc; end eps=abs(b-a)/2; fprintf('\n i=%i a=%f b=%f Ei= %e \n', it, a, b, eps); end x=(a+b)/2; fprintf('\n x=%f \n', x);
Funksione te reja input('f=')
Funksioni i Matlab-it per te lexuar te dhena nga perdoruesi. 14
sign(x)
Funksioni shenje
fprintf('\n i=%i', it);
Funksioni per te printuar ne dritaren e komandave. Nepermjet ketij funksioni realizohet printimi i shume variablave dhe formatimi i variablave. Funksioni nga i cili mund te ndertojme funksione nga stringu.
inline(expression)
Te zbatohet skripti per te gjetur rrenjen e ekuacionit me saktesi .
ne intervalin
Skriptin e mesiperm e shkruajme ne nje skedare me emrin m_pergjysmimit.m dhe nga dritarja e komandave therasim kete skript duke shkruar emrin e tij >> m_pergjysmimit f='x^4+2*x^3-x-1' a=0 b=1 tol=0.01 i=1 a=0.500000 b=1.000000 Ei= 2.500000e-01 i=2 a=0.750000 b=1.000000 Ei= 1.250000e-01 i=3 a=0.750000 b=0.875000 Ei= 6.250000e-02 i=4 a=0.812500 b=0.875000 Ei= 3.125000e-02 i=5 a=0.843750 b=0.875000 Ei= 1.562500e-02 i=6 a=0.859375 b=0.875000 Ei= 7.812500e-03 i=7 a=0.859375 b=0.867188 Ei= 3.906250e-03 x=0.863281
Duke u bazuar ne skriptin e mesiperm te ndertohet funksioni per metoden e pergjysmimit. Te provohet funksioni per ekuacionin e mesiperm per saktesi te ndryshme . Ushtrim. Nepermjet metodes se pergjysmimit perafroni rrenjet e ekuacioneve: a) b)
15
3.2
Metoda e Njutonit
Formula e metodes se Njutoni eshte:
Te ndertohet nje funksion ne Matlab per metoden e Njutonit. Funksioni te marr parasysh rastin e nje funksioni cfaredo si dhe mundesi per te parandaluar rastet kur metoda divergjon. function [ x, eps ] = metoda_njutonit(f, df, x0, tol, nmax) % % Metoda e Njutonit % Metoda e njutonit per gjetjen e rrenjes se ekuacionit % % Hyrjet: % f - funksioni te cilit do ti gjendet rrenja % df - derivati i funksionit % x0 - perafrimi fillestar % tol - toleranca % nmax - numri maksimal i iteracioneve qe i lejojme metodes % % Daljet: % x - rrenja % eps - gabimi ne iteracionin e fundit % % Shembull: % [x, eps] = metoda_njutonit('exp(x)+x', 'exp(x)+1', 0, 0.5*10^-5, 10) if nargin == 3 tol = 1e-4; nmax = 1e1; elseif nargin == 4 nmax = 1e1; elseif nargin ~= 5 error('Metoda e Njutonit: parametra hyres te gabuar!'); end f = inline(f); df = inline(df); eps = tol+1; it = 0; while (eps >= tol) && (it <= nmax) x = x0 - (f(x0)/df(x0)); eps = abs(x-x0); x0=x; it = it+1; fprintf('\n i=%i, x=%f, eps=%e', it, x, eps); end end
nargin
Variabel special qe perdoret vetem ne trupin e funksionit. Ky variabel mban numrin e parametrave qe i kan kaluar funksionit 16
Funksionin e mesiperm mund ta zbatojme per te gjetur rrenjen e ekuacionit , qe ka si derivat , duke marr si perafrim fillestar te rrenjes vleren 0 dhe numer maksimal iteracionesh 10. Ne dritaren e komandave shkruajme >> [ x, eps ] = metoda_njutonit( 'exp(x)+x', 'exp(x)+1', 0, 0.5*10^-5, 10 ) Funksionin metoda_njutonit mund ta therrasim edhe me 4 parametra duke shkruar >> [ x, eps ] = metoda_njutonit( 'exp(x)+x', 'exp(x)+1', 0, 0.5*10^-5 ) Ne kete rast funksioni eshte therritur me 4 parametra. Parametri i fundit eshte 100 si vlere e zbatuar kur nuk percaktohet parametri i fundit. Ne kete menyre ne mund te mos shkruajme parametrin 4 dhe 5 nese duam te njehsojme rrenjen me saktesi dhe numer maksimal iteracionesh 100.
3.3
Metoda e kordes (Lagranzhit)
Ne ngjashmeri me metoden e Njutonit ndertoni nje funksion per metoden e kordes.
3.4
Metoda e pikes fikse
Ndertoni nje funksion per metoden e pikes fikse.
Zbatoni metodat e mesiperme per funksionet e meposhtme per te gjentur rrenjen e tyre dhe krahasoni numrin e iteracioneve te kater metodave. a) b)
17
4. Sisteme te ekuacioneve lineare 4.1
Zgjidhja e sistemeve me matrice trekendeshe te poshtme
Te ndertohet nje algoritem per zgjidhjen e nje sistemi me matrice trekendeshe te poshtme. Formula per zgjidhjen e nje sistemi me matrice trekendeshje e poshtme eshte:
Nga ekuacioni i pare gjendet , duke ditur gjejme duke zevendesuar ne ekuacionin e dyte e keshtu me radhe deri tek ekucaioni i fundit per te gjetur . Kodi ne Matlab ne formen e nje skripti eshte: A=[1, 0, 0; 1, 2, 0; 1, 2, 3]; b=[1, 2, 3]; n=length(A); x=zeros(n, 1); for i=1:n s=0; for j=1:(i-1) s = s + A(i,j)*x(j); end x(i) = (b(i) - s) / A(i,i); end disp(x);
length(X)
Funksioni i gatshem i Matlab-it per te gjetur gjatesin e vektorit kur X eshte vektor ose ne rastin kur X eshte matrice, funksioni na kthen permasen me te madhe te matrices.
Ne ngjashmeri me kodin e mesiperm te ndertohet nje algoritem per zgjidhjen e nje sistemi me matrice trekendeshe te siperme.
4.2
Metoda e Gausit
Te ndertohet nje funksion per metoden e Gausit. Te merret parasysh rasti kur elementi pivot mund te qelloj i barabart me zero. Metoda e Gausit eshte ndertuar ne formen e nje funksioni. Faza e pare e metodes eshte ajo e eliminimit dhe faza e dyte ajo e zevendesimit. Funksioni merr parasysh rastin kur haset element pivon i barabart me zero dhe bene ndrimin me rreshtin e pare qe gjen elementin e pare te ndryshem nga zero nen elementin celes. function [x] = eliminimi_gausit(A,b) %
18
% Zgjidhja e sistemit Ax = b nepermjet eliminimit te Gausit % % Hyrjet: % A - matrica e sistemit % df - vektori i termave te lire % % Daljet: % x - vektori zgjidhje i sistemit Ax=b % gjejme permasat e matrices n = length(b); % faza e eliminimit for j=1:(n-1) if A(j,j) == 0 for i=(j+1):n if A(i,j) ~= 0 tmp = A(j, :); A(j, :) = A(i, :); A(i, :) = tmp; tmp = b(j); b(j) = b(i); b(i) = tmp; break; end end end for i =(j+1):n fprintf('\n i=%i', i); gij = A(i,j)/A(j,j); A(i,j:end) = A(i,j:end)-gij*A(j,j:end) ; b(i) = b(i)-gij*b(j); end end % faza e zevendesimit for i=n:-1:1 s=0; for j=(i+1):n s = s + A(i,j)*x(j); end x(i) = (b(i) - s)/A(i,i); end x = x'; return
Po te ekzekutojme funksionin e siperme per sistemin
Ne dritaren e komandave te Matlab-it shkruajm komandat:
19
ku
>> A=[2,1,0; -1, 4, 2; 1, -1, 3]; >> b=[1, 0, -1]; >> eliminimi_gausit(A,b) ans = 0.3636 0.2727 -0.3636 Pra zgjidhja e sistemi eshte vektori
4.3
.
Faktorizimi LU
Me ndihmen e algoritmit te metodes se Gausit te ndertohet nje funksion qe njehson matricat e faktorizimti .
4.4
Metoda e Jakobit
Te ndretohet nje funksion per zgjidhjen e nje sistemi ekuacionesh lineare me metoden iterative te Jakobit. Formula e metodes se Jakobit eshte:
Ndertome ne Matlab nje funsion qe e emertojme metoda_jakobit. Parametrat hyres dhe dales te funksionit jane shpjeguar ne trupin e funksionit function [x, it]=metoda_jakobit(A,b,x0,tol, nmax) % metoda iterative e Jakobit per zgjidhjen e sistemit Ax=b % hyrjet: % A - eshte matrice n x n permasa % b - eshte vektor me n x 1 permasa % x0 - vektor me n x 1 permasa, permban vetorin fillestar % tol - numer i vogel qe perfaqeson tolerancen % nmax - numri maksimal i iteracioneve % daljet % x - eshte vektor n x 1 permasa, vektori zgjidhje i sistemin Ax=b % it - numri i iteracioneve qe jane kryer per perftimin e x n = length(b); x=zeros(n,1); eps = tol+1; it = 0; while (eps > tol && it < nmax) for i=1:n s = 0;
20
for j=1:n if i ~= j s = s + A(i,j)*x0(j); end end x(i) = (b(i) - s)/A(i,i); end eps = norm(x-x0, inf)/norm(x, inf); it = it +1; x0=x; end end
Nese te dhenat hyrese nuk perputhen sipas pershkrimit tek kodi, atehere mund te shfaqet gabim gjate ekzekutimit te funksionit. Per te mos e komplikuar kodin e mesiperm nuk eshte testuar nese sistemi eshte me diagonale dominante. Megjithate eshte vendosur nje number iteracionesh maksimal qe mund te kryej metoda per te parandaluar divergjencen. Provojme te ekzekutojme kodin e mesiperm per nje sistem ekuacionesh lineare. >> [x, iter]=metoda_jakobit(A, b, zeros(3, 1), 0.0000005, 100) x= 0.186119871221027 0.331230319140125 -0.422712897144088 iter = 13 Zgjidhja e sistemit eshte vektori . Numri i iteracioneve te kryera me metoden e Jakobit per te arritur kete zgjidhje me saktesi 0.0000005 eshte 13 iteracione duke filluar nga vektori fillestar .
Te ndertohet nje funksion per metoden e Gaus-Zajdelit. Te ekzekutohet metoda e Gaus-Zajdelit per sistemin e mesiperm dhe te krahason numri i iteracioneve te kryera midis metodes se Jakobit dhe asaj te Gaus-Zajdelit.
21
5. Interpolimi i funksioneve 5.1
Polinomi interpolues i Lagranzhit
Te ndertohet nje funksion qe njehson vleren e polinomit te Lagranzhit ne nje ose disa pika njeheresh. function [y, L]=Lagranzh(pikat_x, pikat_y, x) % Njehsimi i vleres se polinomit te Lagranzhit ne disa pika % pikat_x, pikat_y jane dy vektor qe perfaqesojne funksionin tabelor % x - vlerat e pikave ne te cilat do gjendet vlera e polinomit % y - vlerat e njehsuara n=size(pikat_x, 2); L=ones(n,size(x, 2)); for i=1:n for j=1:n if (i~=j) L(i,:)=L(i,:).*(x-pikat_x(j))/(pikat_x(i)-pikat_x(j)); end end end y=0; for i=1:n y = y + pikat_y(i)*L(i,:); end end
Ne dritaren e komandave shkruajm >> xi=[0 0.1 0.3 0.5]; >> yi=[-0.5 0 0.2 1]; >> Lagranzh(xi,yi, 0.2) ans = 0.1500 Sic mund te shikoni vlera polinomit interpolues i Lagranzhit ne piken 0.2 eshte 0.15. Kjo vlere eshte ndermjet dy vlerave 0 dhe 0.2 perkatsisht per dhe per . Funksionin e mesiperm mund ta perdorim edhe per te njehsuar disa vlera njeheresh si ne shembullin e meposhtem >> xi=[0 0.1 0.3 0.5]; >> yi=[-0.5 0 0.2 1]; >> Lagranzh(xi,yi, [0.2 0.25]) 22
ans = 0.1500 0.1719 Vlera e polinomit interpolues te lagranzhit ne piken 0.2 eshte 0.15 ndersa ne piken 0.25 eshte 0.1719.
5.2
Polinomi interpolues i Njutonit
Te ndertohet nje funksion qe njehson vleren e polinomit interpolues te Njutonit me diferenca te ndara. Funksioni te mundesoj njehsimin e vleres se polinomit ne disa pika duke perdorur te gjithe informacionin e funksionit tabelor. function [y, D, a] = InterpolimiNjutonit(pikat_x, pikat_y, x) % Njehsimi i vleres se polinomit interpolues i Njutoni me diferenca te % ndara ne nje ose disa pika njeheresh % pikat_x,pikat_y - vektor rresht me gjatesi (n+1) qe perfaqesojne pikat (x,y) ku njihet % vlera e funksioni tabelor % x - vektori me vlerat ku kerkohet te gjendet vlera e polinomit % interpolues % y - vlerat e polinomit interpolues ne pikat x % fuqia e polinomit interpolues n = length(pikat_x) - 1; % numri i vlerave ku kerkohet vlera e polinomit interpolues ni = length(x); D = zeros(n,n); % matrica me diferenca te ndara for k = 1 : n D(k,1) = (pikat_y(k+1)-pikat_y(k))/(pikat_x(k+1)-pikat_x(k)); % diferencat e rendit te pare end for k = 2 : n % cikli per rendin e diferencave te ndara % Njehso te gjitha diferenca te rendit k for i = 1 : (n-k+1) D(i,k) = (D(i+1,k-1) - D(i,k-1))/(pikat_x(i+k)-pikat_x(i)); end end % a(1) = pikat_y(1); for k = 2 : (n+1) a(k) = D(1,k-1); end % Njehsimi i vleres se polinomit interpolues sipas skemes se Hornerit y = a(n+1)*ones(1,ni); for k = 1 : n y = a(n+1-k)+y.*(x-pikat_x(n+1-k)); end end
23
Ekzekutojme ne dritaren e komandave keto komanda >> xi=[0 0.1 0.3 0.5]; >> yi=[-0.5 0 0.2 1]; >> InterpolimiNjutonit(xi, yi, [0.2 0.25]) ans = 0.1500 0.1719 Gjithashtu mund te perftojme edhe diferencat e ndara qe na jep funksioni InterpolimiNjutonit duke shkruar kete komande >> [y, D]=InterpolimiNjutonit(xi, yi, [0.2 0.25]) y= 0.1500 0.1719 D= 5.0000 -13.3333 41.6667 1.0000 7.5000 4.0000
0
0 0
Po te vini re vlerat e poinomit interpolues te Lagranzhit dhe te Njutonit ne pikat 0.2 dhe 0.25 kane dal te njejta. Ketu ka qelluar rastesisht qe jane te njejtat.
24
6. Perafrimi i funksioneve 6.1
Metoda e katroreve me te vegjel
Te ndertohet nje funksion per njehsimin e koeficienteve te polinomit te perafrimit me te mire te rendit kur jepet nje varg pikash . Polinomi i rendit
eshte
Per te gjetur koeficientat e polinomit perafrues me te mire te rendit normale
ku koeficientet
formohet sistemi ekuacioneve
njehsohen sipas formulave
Zgjidhja e sistemit te ekuacioneve lineare me ndonje nga metodat e drejtperdrejta na jep koeficientat e polinomit perafrues te rendit . Ne Matlab ndertojme funksionin qe e quajm PerafrimiPolinomial qe njehson koeficientet e polinomit perafrues te rendit . function [p, A, b]=PerafrimPolinomial(pikat_x, pikat_y, k) % Njehsimi i koeficienteve te polinomit perafrues te rendit k sipas % metodes se katroreve me te vegjel % Hyrjet: % pikat_x, pikat_y - vlerat e dhena te funksionit tabelor % k - rendi i polinomit perafrues % daljet: % p - koeficientet e polinomit perafrues te rendit k % Pk(x)=p(1)+p(2)*x+p(3)*x^2 + ... + p(k+1)*x^k; % A - matrica e sistemit normal % b - vektori i termave te lira te sistemit normal % m - numri i ekuacioneve normale m=k+1; A=zeros(m,m);
25
b=zeros(m,1); % njehsojme elementet e matrices for i=1:m for j=1:i A(i, j)=sum(pikat_x.^((i-1)+(j-1))); A(j, i)=A(i,j); end b(i)=sum(pikat_y.*pikat_x.^(i-1)); end % zgjidhim sistemin e ekuacioneve lineare vektori shtyp p p=A^(-1)*b; % nga vektor shtyll e kthejme ne vektor rresht p=p'; end
Provojme te gjejme polinomin e perafrimit me te mire te rendit te pare per kete funksion tebelor 1 1.3
2 3.5
3 4.2
4 5.0
5 7.0
6 8.8
7 10.1
8 12.5
Ne dritaren e komandave shkruajme keto komanda >> x=1:11; >> y=[1.3 3.5 4.2 5 7 8.8 10.1 12.5 13 15.6 16.1]; >> [p, A, b]= PerafrimPolinomial (x,y, 1) p= -0.2764 1.5173 A= 11 66 66 506 b= 97.1000 749.5000
Polinomi interpolues i rendit te pare eshte
.
Gjithashtu funksioni na jep edhe sistemin e ekuacioneve normale 26
9 13.0
10 15.6
11 16.1
Zgjidhja e ketij sistemi ekuacionesh jep Mund te perftojme edhe polinomin e perafrimit me te mire te rendit te dyte duke shkruar >> PerafrimPolinomial(x,y, 2) ans = 0.0721 1.3564 0.0134 Polinomi interpolues me i mire i rendit te dyte eshte
.
Matlab-i ka nje funksion per perafrimin polinomial qe eshte polyfit. Ky funksion jep koeficientet e polinomit te perafrimit me te mire te rendit k. Shikoni shembullin >> polyfit(x,y, 2) ans = 0.0134 1.3564 0.0721 Pergjigja eshte afishuar duke filluar nga koeficienti i rendit me te lart te polinomit e deri tek koeficienti . Funksioni PerafrimiPolinomial afishon koeficientet ne rendit e kundert nga polyfit.
27
7. Integrimi i perafert i funksioneve 7.1
Metoda e perbere e trapezit
Formula e metodes se perbere te trapezit eshte
Ku
jane pikat e njohura te funksionit qe kerkohet te intergrohet.
Ne Matlab ndertojme funksionin MetodaTrapezit qe kerkon vlerat e funksionit tabelor dhe jep vleren e integralit te njehsuar sipas formules se perbere te trapezit. function [ I ] = MetodaTrapezit(pikat_x, pikat_y) % % Njehsimi i integralit te funksionit tabelor sipas % metodes se perbere te trapezit % % Hyrjet: % pikat_x - vlerat x(i) te baraslagruara % pikat_y - vlerat y(i) ku njihet funksioni % % Daljet: % I - vlera e njehsuar e integralit n=length(pikat_x); a = pikat_x(1); b = pikat_x(n); h=(b-a)/(n-1); I = h * (pikat_y(1) + pikat_y(n))/2 + h*(sum(pikat_y(2:(n-1)))); end
Vini re ndryshimin e indekseve tek funksioni duke krahasuar me ato te formules. Matlab-i indekson vektoret nga 1 ndersa formula e metodes se trapezit indeksi i pare eshte . Provojme te ekzekutojme ne dritaren e komandave >> x=0:0.1:1; >> y=[1, 0.99, 0.9608, 0.9439, 0.8521, 0.7788, 0.6977, 0.6126, 0.5273, 0.4449, 0.3679]; >> MetodaTrapezit(x,y) ans = 0.7492
28
Vlera e integralit te caktuar ne [0, 1] eshte 0.7492.
7.2
Metoda e perbere e Simpsonit
Formula e metodes se perbere te Simpson-it eshte
Ku
jane pikat e njohura te funksionit qe kerkohet te intergrohet.
Ne Matlab ndertojme funksionin qe do e quajm MetodaSimpsonit. function [I, S1, S2, S3] = MetodaSimpsonit(pikat_x, pikat_y) % % Njehsimi i integralit te funksionit tabelor sipas % metodes se perbere te Simpsonit % % Hyrjet: % pikat_x - vlerat x(i) te baraslagruara % pikat_y - vlerat y(i) ku njihet funksioni % % Daljet: % I - vlera e njehsuar e integralit
n=length(pikat_x); if mod(n, 2) == 0 error('Metoda e Simpsonit kerkon nje numer cift segmentesh!'); end a = pikat_x(1); b = pikat_x(n); h=(b-a)/(n-1); % S1 shuma e y(0) + y(n) S1 = pikat_y(1) + pikat_y(n); % S2 shuma e elementeve me index tek S2 = sum(pikat_y(2:2:(n-1))); % S3 shuma e elementeve me index cift S3 = sum(pikat_y(3:2:(n-2))); I = h /3 * (S1 + 4*S2 + 2*S3); end
Duhet te kemi parasysh qe metoda e Simpsonit kerkon nje numer cift segmentesh ose nje numer tek pikash te dhena. Nese nuk plotesohet nje numer cift segmentesh, funksioni do te afishoj mesazhin 'Metoda e Simpsonit kerkon nje numer cift segmentesh!' qe sinjalizon se numri i segmenteve nuk eshte cift.
29
Ne dritaren e komandave shkruajme >> x=0:0.1:0.8; >> y=[1, 0.995, 0.9801, 0.9553, 0.9211, 0.8776, 0.8253, 0.7648, 0.6967]; >> MetodaSimpsonit(x,y) ans = 0.7174 >> MetodaTrapezit(x,y) ans = 0.7168 Vlera qe na jep metoda e perbere e Simpsonit eshte 0.7174. Ndersa ajo e perbere e trapezit na jep vleren 0.7168. Po te duam vlerat e shumes se termave te pare dhe te fundit, tek ose cift mund te percaktojme parametrat dales S1, S2, S3 si meposhte >> [I, s1, s2, s3]=MetodaSimpsonit(x,y) I= 0.7174 s1 = 1.6967 s2 = 3.5927 s3 = 2.7265
30
8. Zgjidhja e perafert e ekuacioneve diferenciale te zakonshme 8.1
Metoda e Eulerit
Te ndertohet nje funksion per zgjidhjen e problemit te vleres fillestare me metoden e Eulerit. Problemi i vleres fillestare eshte
Ku
ndahet ne
segmenteve te barabarta. Te dhenat qe njohim jane
.
Formula e Eulerit eshte
Funksioni ne Matlab eshte function [y, x]=MetodaEulerit(f, a, b, ya, n) % % Metoda e Eulerit per zgjidhje e problemit Koshi % %Hyrjet: % f - funksioni si string 'f' % a, b - skaji i majt dhe ai i djathe % ya - vlera fillestare y(a) % n - [a, b] ndahet ne n segmente te barabarta %Daljet: % % Y is the vector of ordinates f = inline(f); h=(b-a)/n; x=zeros(1,n+1); y=zeros(1,n+1); x=a:h:b; y(1)=ya; for j=1:n y(j+1)=y(j)+h*feval(f,x(j), y(j)); end end
Metoda e Eulerit kerkon te gjitha parametra qe njihen te problemit te vleres fillestare. Pergjigja e kesaj metode eshte vektori me vlerat e njehsuara sipas metodes se Eulerit. Zbatojme funksionin e mesiperm per perafrimin e vleres fillestare
31
Ne dritaren e komandave shkruajme >> [y,x]=MetodaEulerit('-y+x+1', 0,1, 1, 10) y= 1.0000 1.0000 1.0100 1.0290 1.0561 1.0905 1.1314 1.1783 1.2305 1.2874 1.3487 x= 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 Zgjidhja e problemi te vleres fillesatare eshte
.
Ushtrim. Jepet problemi i vleres fillestare me zgjidhje te sakte
Perafroni
me metoden e Eulerit duke marre
.
Te ndertohen rezultatet e perfuara ne te njejtin grafik.
8.2
Metoda Runge-Kutta
Te ndertohet nje funksion per zgjiedhjen e problemit te vleres fillesatare me metoden e Runge-Kutta. Konsiderojme problemi i vleres fillestare
Ku
ndahet ne
segmenteve te barabarta. Te dhenat qe njohim jane
Formula e metodes Runge-Kutta eshte
32
.
Funksioni ne Matlab eshte function [y, x]=MetodaRungeKutta(f, a, b, ya, n) % % Metoda Runge-Kutta per zgjidhje e problemit Koshi % %Hyrjet: % f - funksioni si string 'f' % a, b - skaji i majt dhe ai i djathe % ya - vlera fillestare y(a) % n - [a, b] ndahet ne n segmente te barabarta %Daljet: % Y - vektori me vlerat e perafrta te funksionit y ne pikat x(i) f = inline(f); h=(b-a)/n; x=zeros(1,n+1); y=zeros(1,n+1); x=a:h:b; y(1)=ya; k1=0; k2=0; k3=0; k4=0; for i=1:n k1= h * feval(f, x(i), y(i)); k2= h * feval(f, x(i) + h/2, y(i)+ k1/2); k3= h * feval(f, x(i) + h/2, y(i)+ k2/2); k4= h * feval(f, x(i) + h, y(i)+ k3); y(i+1)=y(i)+(k1 + 2*k2+ 2*k3 + k4)/6; end end
Metoda e Runge-Kutta kerkon te gjitha parametra qe njihen te problemit te vleres fillestare. Pergjigja e kesaj metode eshte vektori me vlerat te peraferta te funksionit ne pikat . Zbatojme metoden Runge-Kutta per perafrimin e vleres fillestare
Ne dritaren e komandave shkruajme >> [y,x]=MetodaRungeKutta('-y+x+1', 0, 1, 1, 10) y= 1.0000 1.0048 1.0187 1.0408 1.0703 1.1065 1.1488 1.1966 1.2493 1.3066 1.3679 33
x= 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 Metoda na ka dhene vektorin
qe eshte vektori me vlerat e funksionit
ne pikat
.
Tani ndertojme grafikun qe permban vlerat e sakta, vlerat e njehsuara sipas metodes se Eulerit dhe vlerat e njehsuara sipas metodes Runge-Kutta nga skripti i meposhtem x=0:0.2:1; y=exp(-x)+x; [y1,x1]=MetodaEulerit('-y+x+1', 0, 1, 1, 10); [y2,x2]=MetodaRungeKutta('-y+x+1', 0, 1, 1, 2); plot(x, y, '-.r', x1, y1, '-.b', x2, y2, '-.g'); leg=legend('y(i)-Vlera Sakte','y(i)-M.E', 'y(i)-MRK4'); set(leg,'Location','NorthWest')
Ky skript do te na afishoj kete grafike
1.4 y(i)-Vlera Sakte y(i)-M.E y(i)-MRK4
1.35 1.3 1.25 1.2 1.15 1.1 1.05 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Grafiku i vlerave te sakta, vlerave te perftuara nga metoda e Eulerit dhe ajo Runge-Kutta per nje problem te vlerave fillestare. Sic mund te shikoni nga grafiku, metoda Runge-Kutta vetem me dy iteracione ka perafruar vlera me te sakta se ajo e Eulerit me 10 iteracione.
34
