Politechnika Lubelska w Lublinie	Laboratorium Metod Numerycznych
		Ćwiczenie nr 9
Imię i Nazwisko: Jakub Machometa	Semestr: III	Grupa: 3.3	Rok akademicki: 2009/2010
Temat: Metody jednokrokowe rozwiązywania równań różniczkowych.		Data wyk.: 16.12.2009r.	Ocena:

Cel ćwiczenia:

Celem ćwiczenia jest zapoznanie z podstawowymi algorytmami rozwiązywania równań różniczkowych metodami Eulera i Rrungego-Kutty RK4

Schemat układu:

Wartości przyjęte do wykonania skryptu:

podaj wartosc napiecia: 5

podaj pojemnosc kondensatora: 0.0004

podaj indukcyjnosc cewki: 0.02

podaj poczatkoe napiecie na kondensatorze uc(0): 1.65

podaj poczatkowy prad w obwodzie i(0): 0.3

Skrypt SciLab:

clear;

clc;

lines(0);

E=input("podaj wartosc napiecia: ");

C=input("podaj pojemnosc kondensatora: ");

L=input("podaj indukcyjnosc cewki: ");

Rk=(2*sqrt(L/C));

uc0=input("podaj poczatkoe napiecie na kondensatorze uc(0): ");

il0=input("podaj poczatkowy prad w obwodzie i(0): ");

wybmet=input(" Wybierz metode..

1 - METODA KLASYCZNA 2 - METODA EULERA I RUNGEGO-KUTTY RK4..

3 - ROZWIAZANIE METODA RK4")

select wybmet;

case 1 then

function uc=udokladnie(t)

uc=E+A1*exp(s1*t)+A2*exp(s2*t)

endfunction

function il=idokladnie(t)

il=C*(A1*s1*exp(s1*t)+A2*s2*exp(s2*t))

endfunction

R=Rk*5;

a=R/(2*L);

wn=1/(sqrt(L*C));

s1=-a-sqrt(a*a-wn*wn);

s2=-a+sqrt(a*a-wn*wn);

A1=(il0-uc0*s2+E*s2)/(s1-s2);

A2=(il0-uc0*s1+E*s1)/(s2-s1);

t=linspace(0,-(5/s1)-5/s2,200);

xset("window",0);

plot2d(t, udokladnie(t),3);

plot2d(t, idokladnie(t),4);

xtitle('przebiegi dla przypadku aperiodycznego','t','uC i',boxed=1)

legends(['uC','i'],[3,4],opt=2)

function w=uakdokladnie(t)

w=E+(A1+(A2*t)).*(exp(-a*t))

endfunction

function z=iakdokladnie(t)

z=-C*(A2-a*A1+A2*t).*exp(-a*t)

endfunction

R=Rk;

a=R/(2*L);

wn=1/(sqrt(L*C));

s=-a;

A1=uc0-E;

A2=il0+a*(uc0-E);

t=linspace(0,6/a,200);

xset("window",1)

plot2d(t, uakdokladnie(t),3);

plot2d(t, iakdokladnie(t),4);

xtitle('przebiegi dla przypadku aperiodycznego krytycznego','t','uC i',boxed=1)

legends(['uC','i'],[3,4],opt=2)

R=Rk/5;

a=R/(2*L);

wn=1/(sqrt(L*C));

s1=-a-sqrt(a*a-wn*wn);

s2=-a+sqrt(a*a-wn*wn);

A1=(il0-uc0*s2+E*s2)/(s1-s2);

A2=(il0-uc0*s1+E*s1)/(s2-s1);

t=linspace(0,5/a,200);

xset("window",2)

plot2d(t, udokladnie(t),3);

plot2d(t, idokladnie(t),4);

xtitle('przebiegi dla przypadku oscylacyjnego','t','uC i',boxed=1)

legends(['uC','i'],[3,4],opt=6)

case 2 then

function w=f1(t,y1,y2)

w=y2;

endfunction

function w=f2(t,y1,y2)

w=-y1/(L*C)-(R*y2)/L+E/(L*C);

endfunction

function [t,y1,y2]=Euler(t0,tk,h,uc0,il0)

N=(tk-t0)/h;

t(1)=t0;

y1(1)=uc0;

y2(1)=il0;

for n=1:N

t(n+1)=t(n)+h;

y1(n+1)=y1(n)+h*f1(t(n),y1(n),y2(n));

y2(n+1)=y2(n)+h*f2(t(n),y1(n),y2(n));

end

endfunction

R=Rk;

tau=(2*L)/R;

[eu_t, eu_y1, eu_y2]=Euler(0, 6*tau, tau/2, uc0, il0);

uC=eu_y1;

i=eu_y2*C;

xset("window",0)

plot2d(eu_t, [uC,i],[3,4] );

xtitle('Metoda Eulera przypadek aperiodyczny krytyczny','t','uC i',boxed=1)

legends(['uC','i'],[3,4],opt=2)

R=Rk*5;

tau=(2*L)/R;

[rk4_t, rk4_y1, rk4_y2]=Euler(0, 6*tau, tau/2, uc0, il0);

uC=eu_y1;

i=eu_y2*C;

xset("window",1)

plot2d(eu_t, [uC,i],[3,4] );

xtitle('Metoda Eulera przypadek aperiodyczny','t','uC i',boxed=1)

legends(['uC','i'],[3,4],opt=2)

R=Rk/5;

tau=(2*L)/R;

[rk4_t, rk4_y1, rk4_y2]=Euler(0, 6*tau, tau/2, uc0, il0);

uC=eu_y1;

i=eu_y2*C;

xset("window",2)

plot2d(eu_t, [uC,i],[3,4] );

xtitle('Metoda Eulera przypadek oscylacyjny','t','uC i',boxed=1)

legends(['uC','i'],[3,4],opt=4)

case 3 then

function w=f1(t,y1,y2)

w=y2;

endfunction

function w=f2(t,y1,y2)

w=-y1/(L*C)-(R*y2)/L+E/(L*C);

endfunction

function [t,y1,y2]=RK4(t0, tk, h, uc0, il0)

N=(tk-t0)/h;

t(1)=t0;

y1(1)=uc0;

y2(1)=il0;

for n=1:N

t(n+1)=t(n)+h;

k11=f1(t(n), y1(n), y2(n));

k21=f2(t(n), y1(n), y2(n));

k12=f1(t(n)+ h/2, y1(n)+h/2*k11, y2(n)+h/2*k21);

k22=f2(t(n)+ h/2, y1(n)+h/2*k11, y2(n)+h/2*k21);

k13=f1(t(n)+ h/2, y1(n)+h/2*k12, y2(n)+h/2*k22);

k23=f2(t(n)+ h/2, y1(n)+h/2*k12, y2(n)+h/2*k22);

k14=f1(t(n)+ h, y1(n)+h*k13, y2(n)+h*k23);

k24=f2(t(n)+ h, y1(n)+h*k13, y2(n)+h*k23);

y1(n+1)=y1(n)+h/6*(k11+2*k12+2*k13+k14);

y2(n+1)=y2(n)+h/6*(k21+2*k22+2*k23+k24);

end

endfunction

R=Rk;

tau=(2*L)/R;

[rk4_t, rk4_y1, rk4_y2]=RK4(0, 6*tau, tau/2, uc0, il0);

uC=rk4_y1;

i=rk4_y2*C;

xset("window",0)

plot2d(rk4_t, [uC,i],[3,4] );

xtitle('Metoda RK4 przypadek aperiodyczny krytyczny','t','uC i',boxed=1)

legends(['uC','i'],[3,4],opt=2)

R=Rk*5;

tau=(2*L)/R;

[rk4_t, rk4_y1, rk4_y2]=RK4(0, 6*tau, tau/2, uc0, il0);

uC=rk4_y1;

i=rk4_y2*C;

xset("window",1)

plot2d(rk4_t, [uC,i],[3,4] );

xtitle('Metoda RK4 przypadek aperiodyczny ','t','uC i',boxed=1)

legends(['uC','i'],[3,4],opt=2)

R=Rk/5;

tau=(2*L)/R;

[rk4_t, rk4_y1, rk4_y2]=RK4(0, 6*tau, tau/2, uc0, il0);

uC=rk4_y1;

i=rk4_y2*C;

xset("window",2)

plot2d(rk4_t, [uC,i],[3,4] );

xtitle('Metoda RK4 przypadek oscylacyjny','t','uC i',boxed=1)

legends(['uC','i'],[3,4],opt=4)

end;

Wykresy wykreślone przez SciLab:

WNIOSKI:

Metoda Eulera jest mniej dokładniejsza od metody Rungego-Kutty rzędu 4, przy niskiej liczbie kroków. Gdy liczba kroków zostanie zwiększona, na wykresach nie widać dużych rozbieżności miedzy obiema metodami.

---