Politechnika Lubelska w Lublinie	Laboratorium Metod Numerycznych
	Ćwiczenie nr 9
Imię i Nazwisko:	Semestr: III
Temat: Metody jednokrokowe rozwiązywania równań różniczkowych.	Data wyk.: 25.01.2014r.

Cel ćwiczenia:

Celem ćwiczenia jest wykorzystanie programu Scilab do obliczania rozwiązań równań różniczkowych zwyczajnych w stanach nieustalonych obwodów elektrycznych.

Schemat układu:

Wartości przyjęte do wykonania skryptu:

Napięcie: 100

Pojemność kondensatora: 0.0001 F

Indukcyjność cewki: 0.8 H

Początkowe napięcie na kondensatorze uc(0): 15

Początkowy prąd w obwodzie i(0): 0.1

Skrypt SciLab:

clear; //czyszczenie pamieci

clc; //czyszczenie konsoli

lines(0);

E=input("podaj wartosc napiecia: "); //pobranie wartości napięcia

C=input("podaj pojemnosc kondensatora: ");//pobranie wartości pojemności

L=input("podaj indukcyjnosc cewki: ");//pobranie wartości indukcyjnsci

Rk=(2*sqrt(L/C)); // obliczanie rezystancji krytycznej

uc0=input("podaj poczatkoe napiecie na kondensatorze uc(0): "); //pobranie wartosci napiecia poczatkowego

il0=input("podaj poczatkowy prad w obwodzie i(0): "); //pobranie wartosci prądu poczatkowego

wybmet=input(" Wybierz metode.. // wybór metody

1 - METODA KLASYCZNA 2 - METODA EULERA..

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;

Wyniki:

Napięcia - metoda RK4

Aperiodyczny kryt.	Aperiodyczny	Oscylacyjny
15. 22.747666 37.545596 52.651035 65.533949 75.606519 83.08689 88.456205 92.218037 94.806353 96.562385 97.74046 98.523584	15. 15.318452 15.966822 16.734725 17.542606 18.360415 19.176834 19.98763 20.791255 21.587173 22.37523 23.15542 23.927788	15. 75.872629 98.915818 102.3643 101.22477 100.2868 99.987416 99.960088 99.982994 99.99686 100.00056 100.00064 100.00023

Metoda Eulera

Aperiodyczny kryt.	Aperiodyczny	Oscylacyjny
15. 15.000447 36.250447 57.500335 73.437724 84.06264 90.703209 94.687549 97.011747 98.339859 99.086923 99.501958 99.730227	15. 15.000447 36.250447 57.500335 73.437724 84.06264 90.703209 94.687549 97.011747 98.339859 99.086923 99.501958 99.730227	15. 15.000447 36.250447 57.500335 73.437724 84.06264 90.703209 94.687549 97.011747 98.339859 99.086923 99.501958 99.730227

Wykresy wykreślone przez SciLab:

Metoda klasyczna:

Metoda Eulera RK4:

Metoda RK4

Wnioski: