Laboratorijske vaje Numeriˇcne metode 7. Vaja. Sistemi linearnih enaˇcb B. Jurˇciˇc Zlobec Numeriˇcne metode FE, Ljubljana, 9. december 2013

Laboratorijske vaje Numeriˇ cne metode

7. Vaja. Sistemi linearnih enaˇ cb

B. Jurˇ ciˇ c Zlobec

Numeriˇ cne metode FE, Ljubljana, 9. december 2013

Robni problem

Poiˇsˇ ci pribliˇ zno reˇsitev robnega problema.

¨

y(t ) + β y(t) = ˙ −g y(0) = 0, y(a) = 0.

interval [0, a] razdeliˇs na n podintervalov.

Krajiˇsˇ ca x 0 = 0, x i , i, 1, . . . , n − 1 in x n = a.

V x _i , i = 1, . . . , n − 1 izraˇ cunamo pribliˇ zno vrednost reˇsitve.

y i−1 − 2y i + y i+1

h ² + β y i+1 − y i−1

2h = −g , kjer je y ₀ = 0, y _n = 0, i = 1, . . . , n − 1.

Vzemi a = 1, n = 15, β = 0.12, g = 9.80665 in h = a/n.

Koliko je najveˇ cja od izraˇ cunanih vrednosti y _i ?

Zapiˇsimo sistem enaˇ cb in ga reˇsimo z levim deljenjem.

Pri tem bomo uporabili ukaz diag.

Octave

n=15;

a=1; b=0.12;

g=9.80665; h=a/n;

l0=1/hˆ2-b/2/h;

d0=-2/hˆ2;

u0=1/hˆ2+b/2/h;

d=ones(n-1,1)*d0;

l=ones(n-2,1)*l0;

u=ones(n-2,1)*u0;

A=diag(d)+diag(l,-1)+diag(u,+1);

B=-ones(n-1,1)*g;

Y=A\B;

max(Y)

Gaussova eliminacija

%resi sistem enacb

A=[2 3 4 5;1 2 1 2;1 3 -2 2; 2 -3 4 5]

b=[2;3;4;1];

x=A\b;

b’

x’

A = 2 3 4 5

1 2 1 2

1 3 -2 2 2 -3 4 5 ans = 2 3 4 1

ans = 9.50000 0.16667 -0.27778 -3.27778

%Gaussova eliminacija C=[A,b]

%

C(1,:)=C(1,:)/C(1,1)

%

C(2,:)=C(2,:)-C(2,1)*C(1,:) C(3,:)=C(3,:)-C(3,1)*C(1,:) C(4,:)=C(4,:)-C(4,1)*C(1,:)

% zamenjamo pivotni vrstici tmp=C(2,:)

C(2,:)=C(4,:) C(4,:)=tmp

%

C(2,:)=C(2,:)/C(2,2)

C(3,:)=C(3,:)-C(3,2)*C(2,:) C(4,:)=C(4,:)-C(4,2)*C(2,:) C(3,:)=C(3,:)/C(3,3)

C(4,:)=C(4,:)-C(4,3)*C(3,:)

C(4,:)=C(4,:)/C(4,4)

1.00000 1.50000 2.00000 2.50000 1.00000 -0.00000 1.00000 -0.00000 -0.00000 0.16667 -0.00000 -0.00000 1.00000 0.12500 -0.68750 -0.00000 -0.00000 -0.00000 1.00000 -3.27778

for j=3:-1:1 for i=j:-1:1

C(i,:)=C(i,:)-C(i,j+1)*C(j+1,:) end;

end;

1.00000 0.00000 0.00000 0.00000 9.50000

0.00000 1.00000 0.00000 0.00000 0.16667

0.00000 0.00000 1.00000 0.00000 -0.27778

-0.00000 -0.00000 -0.00000 1.00000 -3.27778

Program

%%Gaussova-Jordanova eliminacija delno pivotirtanje function [C,p]=gaussova(A,b)

C=[A,b]; n=length(b); p=[];

for i=1:n-1

[a,k]=max(abs(C(i:end,i)));

p=[p,C(k+i-1,i)];

temp=C(k+i-1,:); C(k+i-1,:)=C(i,:); C(i,:)=temp;

C(i,:)=C(i,:)/C(i,i);

for j=i:n-1

C(j+1,:)=C(j+1,:)-C(j+1,i)*C(i,:);

end;

end;

C(n,:)=C(n,:)/C(n,n);

for j=n-1:-1:1 for i=j:-1:1

C(i,:)=C(i,:)-C(i,j+1)*C(j+1,:);

end;

end;

Reˇsi robni problem

y ⁰⁰ (x ) + x y(x ) = 0, y(0) = 0, y(2) = 1,

Pribliˇ zna reˇ sitev

Interval [0, 2] razdelimo na n = 20 enakih delov.

Pribliˇ zne vrednosti funkcije y i = y (x i ) v vozliˇsˇ cih

x i = i h, i = 1, . . . , n − 1, h = 2 n , dobimo kot reˇsitev sistema linearnih enaˇ cb:

y i+1 − 2y i + y i−1

h ² + x i y i = 0

i = 1, . . . , n − 1, y 0 = 0, in y n = 2.

Tridiagonalni sistem My = b

Razˇsirjena matrika sistema: [M|b]















−2 + x ₁ h ² 1 0 · · · 0 −y ₀

1 −2 + x 2 h ² 1 · · · 0 0

. . . .

0 · · · 1 −2 + x n−2 h ² 1 0

0 · · · 0 1 −2 + x n−1 h ² −y _n















Reˇsevanje tridiagonalnega sistema

Diagonala matrike M, d

d _i = −2 + x _i h ² , i = 1, . . . , n − 1, zgornja in spodnja vzporednica l in u,

u = (

n−2 krat

z }| {

1, . . . , 1), l = u.

in desna stran sistema

b = (−y ₀ ,

n−3 krat

z }| { 0, . . . , 0, −y _n ).

Zapiˇsi program, ki bo sprejel ˇstiri vektorje d , u, l, b in bo

podal reˇsitev v vektorju y. Sistem reˇsi s pomoˇ cjo Gaussove

eliminacije ne da bi zapisal polno obliko matrike M sistema.

Algoritem (tridiagonalni Gauss + vzvratno vstavljanje)

Gaussova eliminacija:









d

u

0 · · · 0 b

l

d

u

· · · 0 b

. . . .

0 · · · 0 l

d

b









→









δ

u

0 · · · 0 β

0 δ

u

· · · 0 β

. . . .

0 · · · 0 0 δ

β









δ ₁ = d ₁ , δ _k = d _k − u k−1 l k−1 /δ _k −1 , k = 2, . . . n β ₁ = b ₁ , β _k = b _k − u k−1 l k−1 /δ _k −1 , k = 2, . . . n Vzvratno vstavljanje:

y n = β n /δ n , y _k = (β _k − u _k y _k+1 )/δ _k , k = n − 1, . . . , 1

Algoritem (tridiagonalni Gauss + vzvratno vstavljanje)

Gaussova eliminacija:









d

u

0 · · · 0 b

l

d

u

· · · 0 b

. . . .

0 · · · 0 l

d

b









→









δ

u

0 · · · 0 β

0 δ

u

· · · 0 β

. . . .

0 · · · 0 0 δ

β









δ ₁ = d ₁ , δ _k = d _k − u k−1 l k−1 /δ _k −1 , k = 2, . . . n β ₁ = b ₁ , β _k = b _k − u k−1 l k−1 /δ _k −1 , k = 2, . . . n Vzvratno vstavljanje:

y n = β n /δ n , y _k = (β _k − u _k y _k+1 )/δ _k , k = n − 1, . . . , 1

Algoritem (tridiagonalni Gauss + vzvratno vstavljanje)

Gaussova eliminacija:









d

u

0 · · · 0 b

l

d

u

· · · 0 b

. . . .

0 · · · 0 l

d

b









→









δ

u

0 · · · 0 β

0 δ

u

· · · 0 β

. . . .

0 · · · 0 0 δ

β









δ ₁ = d ₁ , δ _k = d _k − u k−1 l k−1 /δ _k −1 , k = 2, . . . n β ₁ = b ₁ , β _k = b _k − u k−1 l k−1 /δ _k −1 , k = 2, . . . n Vzvratno vstavljanje:

y n = β n /δ n , y _k = (β _k − u _k y _k+1 )/δ _k , k = n − 1, . . . , 1

Program funkcija trisys

% x = trisys(u,d,l,b)

% u : nad, l : pod glavno diagonalo, dolzina je n-1

% d : na glavni diagonali, dolzina je n

% b : svobodni cleni, dolzina je n

% x : resitev sistema, dolzina je n

function x = trisys(u,d,l,b) n = length(b);

for k=2:n,mult=l(k-1)/d(k-1);

d(k)=d(k)-mult*u(k-1);

b(k)=b(k)-mult*b(k-1);

end;

x(n) = b(n)/d(n);

for k=(n-1):-1:1,x(k)=(b(k)-u(k)*x(k+1))/d(k); end;

Program funkcija trisys

% x = trisys(u,d,l,b)

% u : nad, l : pod glavno diagonalo, dolzina je n-1

% d : na glavni diagonali, dolzina je n

% b : svobodni cleni, dolzina je n

% x : resitev sistema, dolzina je n

function x = trisys(u,d,l,b) n = length(b);

for k=2:n,mult=l(k-1)/d(k-1);

d(k)=d(k)-mult*u(k-1);

b(k)=b(k)-mult*b(k-1);

end;

x(n) = b(n)/d(n);

for k=(n-1):-1:1,x(k)=(b(k)-u(k)*x(k+1))/d(k); end;

Program funkcija trisys

% x = trisys(u,d,l,b)

% u : nad, l : pod glavno diagonalo, dolzina je n-1

% d : na glavni diagonali, dolzina je n

% b : svobodni cleni, dolzina je n

% x : resitev sistema, dolzina je n

function x = trisys(u,d,l,b) n = length(b);

for k=2:n,mult=l(k-1)/d(k-1);

d(k)=d(k)-mult*u(k-1);

b(k)=b(k)-mult*b(k-1);

end;

x(n) = b(n)/d(n);

for k=(n-1):-1:1,x(k)=(b(k)-u(k)*x(k+1))/d(k); end;

Gravni program

%priprava podatkov

n=20; x0=0; xend=2; y0=0; yend=1; h=(xend-x0)/n;

d=-2+(1:n-1)*hˆ3;

l=ones(1,n-2); u=l;

b=zeros(size(d));

b(1)=-y0; b(end)=-yend;

%klic trisys y=trisys(u,d,l,b);

%graf

x=linspace(x0,xend,n+1); y=[y0,y,yend];

plot(x,y,’o’);

Gravni program

%priprava podatkov

n=20; x0=0; xend=2; y0=0; yend=1; h=(xend-x0)/n;

d=-2+(1:n-1)*hˆ3;

l=ones(1,n-2); u=l;

b=zeros(size(d));

b(1)=-y0; b(end)=-yend;

%klic trisys y=trisys(u,d,l,b);

%graf

x=linspace(x0,xend,n+1); y=[y0,y,yend];

plot(x,y,’o’);

Slika 1

0.5 1.0 1.5 2.0

0.2 0.4 0.6 0.8 1.0 1.2

Slika : Toˇ cna in pribliˇ zna reˇsitev

Reˇsimo gornji sistem z levim deljenjem z redko matriko

Matriko sistema sestavimo jo ukazom M=sparse(I,J,V)

I=[1:n-1,2:n-1,1:n-2];

J=[1:n-1,1:n-2,2:n-1];

V=[-2+(1:n-1)*hˆ3,ones(1,2*n-4)];

M=sparse(I,J,V);

Reˇsimo sistem z levim deljenjem

b=zeros(n-1,1);

b(1)=-y0;

b(end)=-yend;

Y=M\b;

Y=[y0,Y’,yend];

Reˇsimo gornji sistem z levim deljenjem z redko matriko

Matriko sistema sestavimo jo ukazom M=sparse(I,J,V)

I=[1:n-1,2:n-1,1:n-2];

J=[1:n-1,1:n-2,2:n-1];

V=[-2+(1:n-1)*hˆ3,ones(1,2*n-4)];

M=sparse(I,J,V);

Reˇsimo sistem z levim deljenjem

b=zeros(n-1,1);

b(1)=-y0;

b(end)=-yend;

Y=M\b;

Y=[y0,Y’,yend];

Reˇsimo gornji sistem z levim deljenjem z redko matriko

Matriko sistema sestavimo jo ukazom M=sparse(I,J,V)

I=[1:n-1,2:n-1,1:n-2];

J=[1:n-1,1:n-2,2:n-1];

V=[-2+(1:n-1)*hˆ3,ones(1,2*n-4)];

M=sparse(I,J,V);

Reˇsimo sistem z levim deljenjem

b=zeros(n-1,1);

b(1)=-y0;

b(end)=-yend;

Y=M\b;

Y=[y0,Y’,yend];