LABOLATORIUM METOD NUMERYCZNYCH

Ćwiczenie Nr 9: Metody jednokrokowe rozwiązywania równań różniczkowych

Imię i nazwisko: Mateusz Mielniczuk

Mateusz Litwin Gr. 3.3

Wyprowadzenie równań dla obwodu szeregowego RLC:

Rozróżniamy tu trzy przypadki:

1)
_, czyli
_- przypadek aperiodyczny:

,gdzie

2)
_, czyli
_- przypadek aperiodyczny krytyczny:

,gdzie

3)
_, czyli
_- przypadek oscylacyjny:

Wzory zastosowane do obliczenia w scilabie są takie same jak dla przypadku aperiodycznego

METODA KLASYCZNA

clear;

clc;

lines(0);

//parametry zadania

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): ");

//funkcje opisujące przebieg napiecia na kondensatorze i prad w obwodzie

//dla przypadku aperiodycznego i oscylacyjnego

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

//wzory dla przypadku aperiodycznego

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);

xdel(winsid()); //zamkniecie starych wykresów

//wykres dla przypadku aperiodycznego

subplot(221)

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)

//funkcje opisujące przebieg napiecia na kondensatorze i prad w obwodzie

//dla przypadku aperiodycznego krytycznego

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

//wzory dla przypadku aperiodycznego krytycznego

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);

//wykres dla przypadku aperiodycznego krytycznego

subplot(222)

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)

//wzory dla przypadku osylacyjnego

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);

// i wykres

subplot(223)

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)

METODA EULERA I RUNGEGO-KUTTY RZ. 4

Skrypt wpisany w scilabie:

clear; xdel; clc;

lines(0) // wyłączenie przewijania pionowego

//parametry zadania

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

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

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

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

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

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

//rownanie II rzedu zapisane jako uklad rownan rownan I rzedu

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

// ROZWIAZANIE METODA EULERA

//funkcja dla metody Eulera

function [t,y1,y2]=Euler(t0,tk,h,y10,y20)

N=(tk-t0)/h; // obliczanie liczby kroków

t(1)=t0; // inicjacja wektora czasu

y1(1)=y10; // inicjacja wektora rozwiązań

y2(1)=y20; // inicjacja wektora rozwiązań

for n=1:N

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

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

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

end // koniec pętli for,

endfunction // koniec definicji funkcji.

//rozwiązanie dla przypadku aperiodycznego krytycznego

R=Rk;

tau=(2*L)/R;

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

uC=eu_y1; // obliczanie napięcia na kondensatorze

i=eu_y2*C; // obliczanie prądu płynącego w obwodzie

//wykres dla przypadku aperiodycznego krytycznego

subplot(232); // Definicja położenia wykresu

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

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

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

//rozwiązanie dla przypadku aperiodycznego

R=Rk*5;

tau=(2*L)/R;

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

uC=eu_y1; // obliczanie napięcia na kondensatorze

i=eu_y2*C; // obliczanie prądu płynącego w obwodzie

//wykres dla przypadku aperiodycznego

subplot(231); // Definicja położenia wykresu

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

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

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

//rozwiązanie dla przypadku oscylacyjnego

R=Rk/5;

tau=(2*L)/R;

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

uC=eu_y1; // obliczanie napięcia na kondensatorze

i=eu_y2*C; // obliczanie prądu płynącego w obwodzie

//wykres dla przypadku oscylacyjnego

subplot(233); // Definicja położenia wykresu

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

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

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

//ROZWIAZANIE METODA RK4

function [t,y1,y2]=RK4(t0, tk, h, y10, y20)

N=(tk-t0)/h; // obliczanie liczby kroków

t(1)=t0; // inicjacja wektora czasu

y1(1)=y10; // inicjacja wektora rozwiązań

y2(1)=y20; // inicjacja wektora rozwiązań

for n=1:N // rozpoczęcie pętli for,

t(n+1)=t(n)+h; // obliczanie liczby kroków,

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

//rozwiązanie dla przypadku aperiodycznego krytycznego

R=Rk;

tau=(2*L)/R;

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

uC=rk4_y1; // obliczanie napięcia na kondensatorze,

i=rk4_y2*C; // obliczanie prądu w obwodzie,

subplot(235); // Definicja położenia wykresy,

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

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

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

//rozwiązanie dla przypadku aperiodycznego

R=Rk*5;

tau=(2*L)/R;

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

uC=rk4_y1; // obliczanie napięcia na kondensatorze,

i=rk4_y2*C; // obliczanie prądu w obwodzie,

subplot(234); // Definicja położenia wykresy,

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

xtitle('Metoda RK4 przypadek aperiodyczny ','t','uC i',boxed=1) // podpis wykresu.

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

//rozwiązanie dla przypadku oscylacyjnego

R=Rk/5;

tau=(2*L)/R;

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

uC=rk4_y1; // obliczanie napięcia na kondensatorze,

i=rk4_y2*C; // obliczanie prądu w obwodzie,

subplot(236); // Definicja położenia wykresy,

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

xtitle('Metoda RK4 przypadek oscylacyjny','t','uC i',boxed=1) // podpis wykresu.

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

Wykresy dla danych:

E=5V L=0.6 H C=0.05F u_C(0)=0 i_L(0)=0

Wykresy metody klasycznej

Wykresy dla metod: Eulera oraz RK4

WNIOSKI:

Metoda Eulera jest mniej dokładna 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.

---