111
>> pd = polyder(p)
pd =
-19.6217 20.1293
Vi ph©n cña ®a thøcy=-9.8108x2+20.1293x-0.0317lµdx/dy= -19.6217x+20.1293.
Bëi v× ®¹o hµm cña mét ®a thøc còng ®−îc vÏ vµ tÝnh gi¸ trÞ gièng nh− lµ ®èi víi ®a thøc:
>> z = polyval(pd,xi); % evaluate derivative
>> plot(xi,z)
>> xlabel('x'),ylabel('dy/dx')
>> title('Derivative of a Curve Fit Polynomial')
H×nh 16.6
H×nh 16.7
Trong tr−êng hîp nµy xÊp xØ ®a thøc lµ mét hµm bËc hai vµ ®¹o hµm cña nã trë thµnh hµm bËc
nhÊt.
Simpo PDF Merge and Split Unregistered Version - http://www.simpopdf.com
112
MATLAB cung cÊp mét hµm ®Ó tÝnh to¸n ®¹o hµm mét c¸ch s¬ bé dùa vµo d÷ liÖu m« t¶ mét sè
hµm, hµm nµy cã tªn lµ
diff
, nã tÝnh to¸n ®é chªnh lÖch gi÷a c¸c phÇn tö trong m¶ng. Bëi v× ®¹o hµm
®−îc ®Þnh nghÜa nh− sau:
nªn ®¹o hµm cña hµm f(x) cã thÓ ®−îc tÝnh mét c¸ch s¬ bé dùa vµo c«ng thøc:
khi h>0
Gäi lµ sè ra cña y chia cho sè ra cña x, do hµm
diff
tÝnh to¸n sù kh¸c nhau gi÷a c¸c phÇn tö trong
m¶ng nªn ®¹o hµm cã thÓ ®−îc tÝnh mét c¸ch xÊp xØ dùa vµo hµm
diff
:
>> dy = diff(y)./diff(x);
>> % compute differences and use array division
>> xd = x(1:length(x)-1);
>> % create new x axis array since dy is shorter than y
>> plot(xd,dy)
>> title('Approximate Derivative Using DIFF')
>> ylabel('dy/dx'),xlabel('x')
H×nh 16.8
Do hµm
diff
tÝnh ra sù kh¸c nhau gi÷a c¸c phÇn tö nªn kÕt qu¶ cña vÝ dô trªn lµ mét m¶ng cã sè
phÇn tö Ýt h¬n m¶ng ban ®Çu mét phÇn tö. V× vËy ®Ó vÏ ®−îc ®å thÞ cña ®¹o hµm th× ph¶i bá ®i mét
phÇn tö cña m¶ng x. So s¸ng hai ®å thÞ cuèi cïng th× thÊy hiÓn nhiªn r»ng ®¹o hµm tÝnh b»ng ph¬ng
ph¸p gÇn ®óng kh¸c xa so víi thùc tÕ.
16.6 Ph−¬ng tr×nh vi ph©n
Cã thÓ b¹n ®· kh¸ quen víi thùc tÕ lµ rÊt nhiÒu hÖ thèng vËt lý ®Òu ®−îc m« t¶ b»ng ph−¬ng tr×nh
vi ph©n. Do vËy phÇn sau ®©y ®èi víi b¹n cã thÓ kh¸ hÊp dÉn.
Mét ph−¬ng tr×nh vi ph©n th−êng m« t¶ tèc ®é thay ®æi cña mét biÕn sè trong hÖ thèng theo sù
thay ®æi cña mét biÕn kh¸c trong hÖ thèng hoÆc theo kÝch thÝch bªn ngoµi. Ph−¬ng tr×nh vi ph©n th«ng
th−êng cã thÓ ®−îc gi¶i nhê c¸c ph−¬ng ph¸p gi¶i tÝch hoÆc sö dông c«ng cô to¸n kÝ hiÖu cña
MATLAB.
Simpo PDF Merge and Split Unregistered Version - http://www.simpopdf.com
113
Trong nh÷ng tr−êng hîp mµ ph−¬ng tr×nh vi ph©n kh«ng thÓ gi¶i ®−îc b»ng ph−¬ng ph¸p gi¶i tÝch
th× viÖc sö dông ph−¬ng ph¸p sè häc trë lªn kh¸ hiÖu qu¶. §Ó minh ho¹ h·y xÐt ph−¬ng tr×nh Van Der
Pol, ph−¬ng tr×nh biÓu diÔn mét bé dao ®éng.
TÊt c¶ c¸c ph−¬ng ph¸p to¸n häc ®Ó gi¶i ph−¬ng tr×nh d¹ng nµy ®Òu sö dông mét ph−¬ng tr×nh vi
ph©n cao cÊp h¬n, t−¬ng ®−¬ng víi mét tËp ph−¬ng tr×nh vi ph©n bËc nhÊt. §èi víi ph−¬ng tr×nh vi
ph©n trªn th× c¸ch gi¶i nµy ®−îc thùc hiÖn b»ng c¸ch ®Þnh nghÜa hai biÕn trung gian:
®Æt y1 = x, vµ y2 =
suy ra:
§èi víi c¸c hÖ ph−¬ng tr×nh nh− thÕ nµy MATLAB cung cÊp mét tËp c¸c hµm ODE ®Ó gi¶i xÊp xØ
ho¸ chóng mét c¸ch sè häc. Trong quyÓn h−íng dÉn nµy chóng ta kh«ng cã kh¶ n¨ng ®Ó nªu hÕt
nh÷ng néi dung vµ øng dông cña tõng hµm trong bé ODE. §Ó t×m hiÓu thªm vÒ c¸c hµmm ODE øng
dông trong rÊt nhiÒu bµi to¸n thÝ dô, h·y gâ >> odedemo t¹i dÊu nh¾c cña MATLAB. Tr−íc hÕt
chóng ta h·y xÐt vÝ dô sau ®©y, chÝnh lµ vÝ dô
ode45
. Chóng ta ph¶i viÕt mét hµm M_file tr¶ vÒ c¸c
®¹o hµm nÕu biÕt tr−íc c¸c gi¸ trÞ tøc thêi cña y1 vµ y2. Trong MATLAB c¸c ®¹o hµm ®−îc cho bëi
c¸c vector cét, trong tr−êng hîp nµy gäi lµ yprime. T−¬ng tù y1 vµ y2 ®îc viÕt d−íi d¹ng vector cét
y. KÕt qu¶ cña mét hµm M_file nh− sau:
function yprime=vdpol(t,y);
% VDPOL(t,y) returns the state derivatives of
% the Van der Pol equation:
%
% x''-mu*(1-x^2)*x+x=0
%
% let y(1)=x and y(2)=x'
%
% then y(1)'=y(2)
% y(2)'=mu*(1-y(1)^2)*y(2)-y(1)
mu=2; % choose 0< mu < 10
yprime=[y(2)
mu*(1-y(1)^2)*y(2)-y(1)]; % output must be a column
Gi¶ sö thêi gian kÐo dµi tõ 0 ®Õn 30 gi©y, vÝ dô tspan=[0 30]. Sau ®ã sö dông lÖnh
vdpol
th×
lêi gi¶i cho bµi to¸n nh− sau:
>> tspan = [0 30];
>> yo = [1;0];
>> ode45('vdpol',tspan,yo);
Khi sö dông hµm mµ kh«ng cã ®èi sè ra, c¸c hµm ODE sÏ tù ®éng chän nh÷ng thêi ®iÓm thÝch hîp ®Ó
tÝnh ®¹o hµm. §Ó cã thÓ truy nhËp ®−îc d÷ liÖu, ta chØ cÇn cung cÊp cho hµm nh÷ng th«ng sè ra.
>> [t,y] = ode45('vdpol',tspan,yo);
ë ®©y t lµ mét vector cét chøa nh÷ng thêi ®iÓm ®Ó tÝnh ®¹o hµm, cßn y lµ mét ma trËn chøa hai cét
vµ c¸c hµng length(t), hµng ®Çu tiªn cña ma trËn y chøa biÕn sè y(1), hµng thø hai lµ biÕn sè
y(2).
Dùa vµo nh÷ng ®Æc ®iÓm nµy chóng ta cã thÓ vÏ ®−îc ®å thÞ pha, lµ ®å thÞ gi÷a y(2) vµ y(1):
>> plot(y(:,1),y(:,2))
0 5 10 15 20 25 30
-4
-3
-2
-1
0
1
2
3
4
H×nh 16.9
C¸c hµm ODE cña MATLAB ®Òu cã trî gióp trùc tuyÕn, mçi hµm ®Òu cã c¸c ®èi sè còng nh−
c¸ch sö dông riªng, nÕu b¹n muèn nghiªn cøu thªm th× h·y tham kh¶o thªm phÇn trî gióp trùc tuyÕn
cña chóng.
-3 -2 -1 0 1 2 3
-4
-3
-2
-1
0
1
2
3
4
H×nh 16.10
Simpo PDF Merge and Split Unregistered Version - http://www.simpopdf.com
114
ch−¬ng 17
®å ho¹ trong hÖ to¹ ®é ph¼ng
Trong toµn bé phÇn h−íng dÉn sö dông cña cuèn s¸ch nµy, mét sè ®Æc tÝnh vÒ ®å ho¹ cña
MATLAB sÏ lÇn l−ît ®−îc giíi thiÖu, vµ trong ch−¬ng nµy vµ ch−¬ng tiÕp theo chóng ta sÏ lµm s¸ng
tá thªm vÒ nh÷ng ®Æc tÝnh ®ã cña MATLAB.
17.1 Sö dông lÖnh Plot
Nh− b¹n ®· thÊy ë vÝ dô tr−íc ®ã, phÇn lín c¸c c©u lÖnh ®Ó vÏ ®å thÞ trong mÆt ph¼ng ®Òu lµ lÖnh
plot
.LÖnh
plot
nµy sÏ vÏ ®å thÞ cña mét m¶ng d÷ liÖu trong mét hÖ trôc thÝch hîp, vµ nèi c¸c ®iÓm
b»ng ®−êng th¼ng. D−íi ®©y lµ mét vÝ dô mµ b¹n ®· thÊy tr−íc ®ã (H×nh 17.1):
>> x = linspace(0,2*pi,30);
>> y = sin(x);
>> plot(x,y)
VÝ dô nµy t¹o 30 ®iÓm d÷ liÖu trong ®o¹n 0 x 2 theo chiÒu ngang ®å thÞ, vµ t¹o mét vector y kh¸c
lµ hµm sine cña d÷ liÖu chøa trong x. LÖnh
plot
më ra mét cöa sæ ®å ho¹ gäi lµ cöa sæ figure, trong
cöa sæ nµy nã sÏ t¹o ®é chia phï hîp víi d÷ liÖu, vÏ ®å thÞ qua c¸c ®iÓm, vµ ®å thÞ ®−îc t¹o thµnh bëi
viÖc nèi c¸c ®iÓm nµy b»ng ®−êng nÐt liÒn. C¸c thang chia sè vµ dÊu ®−îc tù ®éng cËp nhËt vµo, nÕu
nh−
cöa sæ figure ®· tån t¹i,
plot
xo¸ cöa sæ hiÖn thêi vµ thay vµo ®ã lµ cöa sæ míi.
H×nh 17.1
B©y giê cïng vÏ hµm sine vµ cosine trªn cïng mét ®å thÞ
>> z = cos(x);
>> plot(x,y,x,z)
Simpo PDF Merge and Split Unregistered Version - http://www.simpopdf.com
115
H×nh 17.2
VÝ dô nµy cho thÊy b¹n cã thÓ vÏ nhiÒu h¬n mét ®å thÞ trªn cïng mét h×nh vÏ, b¹n chØ viÖc ®a
thªm vµo
plot
mét cÆp ®èi sè,
plot
tù ®éng vÏ ®å thÞ thø hai b»ng mµu kh¸c trªn mµn h×nh. NhiÒu ®-
−êng cong cã thÓ cïng vÏ mét lóc nÕu nh− b¹n cung cÊp ®ñ c¸c cÆp ®èi sè cho lÖnh
plot
.
NÕu nh− mét trong c¸c ®èi sè lµ ma trËn vµ ®èi sè cßn l¹i lµ vector, th× lÖnh
plot
sÏ vÏ t−¬ng øng
mçi cét cña ma trËn víi vector ®ã:
>> W = [y;z] % x©y dùng mét ma trËn sine vµ cosine
>> plot(x,W) % vÏ c¸c cét cña W víi x
H×nh 17.3
Simpo PDF Merge and Split Unregistered Version - http://www.simpopdf.com