Ćwiczenie Nr 9: Metody jednokrokowe rozwiązywania równań różniczkowych
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 uC(0)=0 iL(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.
