Ćwiczenie 5. Matlab  równania różniczkowe, aproksymacja
Równania różniczkowe - funkcja dsolve()
Funkcja oblicza symbolicznie rozwiązania równań różniczkowych zwyczajnych. Równania są określane przez
symboliczne wyrażenia zawierające literę D do oznaczenia stopnia. Symbole D2, D3... DN, odnoszą się do drugiej,
trzeciej,..., n-tej pochodnej. D2y jest zatem odpowiednikiem symbolicznym . Zmienna niezależna domyślna
to t.
Nazwy zmiennych symbolicznych nie powinny zatem zawierać litery D. Zmienną niezależną można zmienić
i podać jako ostatni argument. Warunki początkowe mogą być określone przez dodatkowe równania. Jeśli nie
określono warunków początkowych, rozwiązania zawierają stałe całkowania: C1, C2, itp.
Przykład 1. Rozwiązać równanie różniczkowe:
( )
= + ( )
f=dsolve('Dy=1+y^2')
=dsolve
=dsolve
=dsolve
f =
tan(t+C1)
f=dsolve('Dy=1+y^2','x')
=dsolve
=dsolve
=dsolve
f =
tan(x+C1)
Po wstawieniu przykładowego warunku początkowego:
f = dsolve('Dy=1+y^2','y(0)=1')
dsolve
dsolve
dsolve
%sprawdzenie
spr= diff(f,t)-1-f^2
f =
tan(t+1/4*pi)
s=
0
Przykład 2. Rozwiązać równanie:
( )
= - ( )
r=dsolve
r=dsolve
r=dsolve('Dx = -a*x')
r=dsolve
r =
C1*exp(-a*t)
Przykład 3. Równanie różniczkowe drugiego stopnia:
( )
( )
= -
z dwoma warunkami poczÄ…tkowymi: y(0)=1; y'(0)=0
y = dsolve('D2y=cos(2*x)-y','y(0)=1','Dy(0)=0', 'x')
dsolve
dsolve
dsolve
y =
4/3*cos(x)-1/3*cos(2*x)
simplify(y) %uproszczenie
simplify
simplify
simplify
ans =
4/3*cos(x)-2/3*cos(x)^2+1/3
Przykład 4. Rozwiązać równanie:
= warunki poczÄ…tkowe: u(0)=1, u'(0)=-1, u"(0)=Ä„
1
u = dsolve
dsolve('D3u=u','u(0)=1','Du(0)=-1','D2u(0) = pi','x')
dsolve
dsolve
u =
1/3*pi*exp(x)-1/3*(pi+1)*3^(1/2)*exp(-1/2*x)*sin(1/2*3^(1/2)*x)+
(1-1/3*pi)*exp(-1/2*x)*cos(1/2*3^(1/2)*x)
D3u reprezentuje d3u/dx3 , D2u(0) odpowiada u"(0)
Układ równań różniczkowych
Funkcja dsolve rozwiązuje także układ równań różniczkowych zwyczajnych kilku zmiennych, z warunkami po-
czÄ…tkowymi lub bez.
Przykład
Dwa równania liniowe, pierwszego rzędu:
( )
( ) ( ) ( ) ( )
= 3 + 4 ; = -4 + 3
S = dsolve
dsolve('Df = 3*f+4*g', 'Dg = -4*f+3*g')
dsolve
dsolve
S =
f: [1x1 sym]
g: [1x1 sym]
Rozwiązania są zwracane w strukturze S. Można określić wartości f i g, wpisując:
f = S.f
f =
exp(3*t)*(C1*sin(4*t)+C2*cos(4*t))
g = S.g
g =
exp(3*t)*(C1*cos(4*t)-C2*sin(4*t))
Jeśli chcemy uzyskać f i g bezpośrednio, oraz uwzględnić także warunki początkowe, wpisujemy:
[f,g] = dsolve('Df=3*f+4*g, Dg =-4*f+3*g', 'f(0) = 0, g(0) = 1')
dsolve
dsolve
dsolve
f =
exp(3*t)*sin(4*t)
g =
exp(3*t)*cos(4*t)
Poniżej jeszcze jeden przykład składni w Symbolic Math Toolbox.
( )
+ 4 = ; y(0)=1
y = dsolve('Dy+4*y = exp(-t)', 'y(0) =1')
dsolve
dsolve
dsolve
spr=diff
diff(y,t)+4*y %sprawdzenie rozwiÄ…zania
diff
diff
spr = simplify
simplify(spr)
simplify
simplify
t=0; %sprawdzenie warunku poczÄ…tkowego
wp=subs
subs(y)
subs
subs
y =
(1/3*exp(3*t)+2/3)*exp(-4*t)
spr =
exp(-4*t)*exp(3*t)
spr =
exp(-t)
wp =
1
2
Zadanie
Równanie oscylatora harmonicznego bez tłumienia to:
( )
+ É0 = 0
gdzie: y(t)  poÅ‚ożenie ciaÅ‚a, É  czÄ™stość drgaÅ„,
zaś oscylatora z tłumieniem:
( )
+ 2² + É = 0
gdzie: ² współczynnik tÅ‚umienia.
Rozwiązać obydwa równania przyjmując różne warunki początkowe. Utworzyć wykresy y(t).
Aproksymacja
Metoda aproksymacji polega na znalezieniu funkcji f(x), której wykres przechodzi w pobliżu zbioru zadanych
punktów, zaś sama funkcja minimalizuje wartość pewnego kryterium dopasowania, najczęściej sformułowanego
w postaci:
[ ( )]
= -
gdzie to wartości zmiennej niezależnej, to wartości zmiennej zależnej. Tego typu aproksymacja nazywa-
na jest aproksymacją średniokwadratową.
Niech w naszym przypadku funkcją aproksymującą f(x) będzie wielomian n-tego stopnia w(x). Wtedy do znale-
zienia współczynników wielomianu używamy funkcji polyfit, której składnia ma postać:
a=polyfit
polyfit(x,y,n)
polyfit
polyfit
gdzie:
x - jest zbiorem N posortowanych rosnąco wartości współrzędnych zmiennej niezależnej,
y - jest zbiorem N odpowiadających wartości zmiennej zależnej,
n - jest zadanym stopniem wielomianu aproksymujÄ…cego,
a - jest wektorem poszukiwanych wartości współczynników wielomianu aproksymującego.
Ćwiczenie
Obliczyć współczynniki paraboli y=p1x2+p2x+p3, która przyjmuje wartości 0 na końcach przedziału <0,10>, a
wartość 1 w połowie przedziału.
Z warunków zadania wynika, że:
x =0, y =0
1 1
x =5, y =1
2 2
x =10, y =0
3 3
Zapisujemy wartości współrzędnych x1,x2,x3 wektorze x, a wartości y1,y2,y3 w wektorze y
x=[0, 5, 10]
y=[0,1,0]
Wykonujemy funkcjÄ™ polyfit
p=polyfit(x,y,2)
Metoda obliczeń w Matlab-ie:
% Obliczenie współczynników paraboli y(x)=p(1)*x^2 + p(2)*x+p(3)
x=[0,5,10];
y=[0,1,0];
p=polyfit(x,y,2)
plot(x,y,'o') %wykres punktowy
% sprawdzenie na wykresie
x1=0:0.1:10;
y1=p(1)*x1.^2+p(2)*x1+p(3);
3
hold on
plot(x1,y1)
Ćwiczenie
Obliczyć współczynniki wielomianu 5-go stopnia, stanowiącego aproksymację dla zbioru punktów:
p1(0, 0), p2(1, 1) , p3(2, 1) , p4(3, 1) , p5(4, 1), p6(5,0)
Przeanalizować poniższe rozwiązanie:
clear
clc
N=5 %stopień wielomianu
x1=[0 1 2 3 4 5]
y1=[0 1 1 1 1 0]
p=polyfit
polyfit(x1,y1,N)
polyfit
polyfit
x=0:0.1:5;
y=0;
%pętla sumująca elementy wielomianu
for m=1:N+1
y=y+p(m)*x.^(N-m+1);
end;
%dwa wykresy
plot(x1,y1,'o',x,y)
axis([0 5 -1 2])
Obliczenie wartości wielomianu o znanych współczynnikach
Wartości otrzymanego wielomianu, w dowolnych zadanych punktach wektorem xx oblicza funkcja polyval.
ypi = polyval
polyval(p,xx)
polyval
polyval
gdzie
p - wektor współczynników wielomianu
ypi - wektor wartości wielomianu w zadanych punktach xx
Realizacja obliczeń:
xx=[0:0.1:10];
p= [1 3 1];
ypi=polyval(p,xx);
plot(xx,ypi) % rysowanie wykresu paraboli w przedziale (0,10)
Ćwiczenie
Zbadać wielomiany stopnia od n=5 do 7, aproksymujące losowo wybrane wartości y z przedziału (0,10) dla
punktów x=[1 2 3 4 5] oraz zgromadzić wykresy funkcji wielomianowych w jednym układzie współrzęd-
nych.
Realizacja obliczeń:
clc
x=1:5;
y=round(10*rand(1,5))
plot(x, y, 'o')
for k=4:7 % pętla stopni wielomianu
hold on
p=polyfit
polyfit(x,y,k)
polyfit
polyfit
xx=0:0.1:5;
yy=polyval
polyval(p,xx);
polyval
polyval
4
 plot(xx,yy)
end
Jak można zauważyć jednoznaczność rozwiązania istnieje dla stopnia wielomianu n mniejszego od liczby punk-
tów N.
Zadania
1. Napisać program, który wczytuje stopień wielomianu n instrukcją:
n = input('n = ')
oblicza współczynniki wielomianu instrukcją polyfit dla danych:
x=[1:12];
y=[1.2 3.7 4.7 4.5 3.4 1.9 1.7 1.8 2.3 3.8 5.4 7.5];
oblicza wartości wielomianu instrukcją polyval dla danych:
xx =[1:0.2:12];
oraz wyświetla wykres instrukcją plot.
2. Obliczyć współczynniki paraboli z=q1y2+q2y+q3, która przyjmuje wartości 0 na końcach przedziału
<0,20>, a wartość 1 w połowie przedziału i narysować wykres paraboli.
5