13 Další možnosti MATLABu a SIMULINKu

13 Další možnosti MATLABu a SIMULINKu Studijní cíl Poslední blok je věnován některým dalším možnostem využití MATLABu a SI‐ MULINKu. V první části si ukážeme možnost volat SIMULINKový model z MATLABu. Druhá část je věnována problematice symbolických výpočtů a vý‐ počtů s požadovanou přesností v MATLABu. Tyto funkce a příkazy nejsou sou‐ částí základní instalace MATLABu, ale jsou obsaženy v rozšiřujícím toolboxu s názvem Symbolic Toolbox. Budeme se stručně zabývat některými příkazy pro řešení úloh z některých oblastí matematiky – lineární algebry, algebraických rovnic, derivací a integrálů funkcí, integrálních transformací a diferenciálních rovnic.

Doba nutná k nastudování

1 hodina

Průvodce studiem Předpokládá se, že čtenář je seznámen s prostředím a možnostmi MATLABu i SIMULINKu popsanými v předchozích kapitolách. Pro řešení příkladů jsou též potřeba základní znalosti matematiky a z oblasti diferenciálních rovnic. Při stu‐ diu je vhodné mít spuštěný MATLAB a SIMULINK a jednotlivé příkazy či funkce hned vyzkoušet a případně získat podrobnější informace využíváním nápovědy k programu. Použití funkcí z Symbolic Math Tbx vyžaduje mít tento Tbx nain‐ stalovaný (není součástí standardní instalace MATLABu). V tomto bloku nejsou uvedeny příklady na procvičení. V textu jsou použity následující typografické konvence: Calibri 11, Bold, Italic Calibri 11 Courier Courier Courier Courier

New New New New

Courier New

nové pojmy k zapamatování označení klávesy 10, Bold názvy nástrojů MATLABu 10 názvy bloků použitých v programu 10 označení části příkazu, která se může vynechat 10, Italic názvy proměnných použitých v programu názvy parametrů bloků či dialogových oken 10, Bold upřesnění nápovědy (help téma)

13.1 Volání modelu SIMULINKu z MATLABu Spolupráce mezi SIMULINKem a MATLABEM nespočívá jen v používání příkazů MATLABu v modelu SIMULINKu a předávání parametrů a průběhů vstupních a výstupních signálů jak bylo popsáno dříve, ale je možné volat SIMULINKový model z příkazové řádky MATLABu. Je možné (bez nutnosti specifikovat v modelu) zadat nejen průběhy vstupních signálů a určit proměnné, do kterých KŘP/IMSW Modelování ve výpočtových software

13‐1 (17) 17.12.11

František Dušek KŘP FEI Univerzita Pardubice

se mají uložit průběhy výstupních signálů, ale je možné zadat počáteční hodno‐ ty a získat hodnoty koncových stavů příslušných bloků (především bloku Integrator). Je také možné nastavit všechny parametry řešení a spuštění modelu včetně doby řešení. Volání SIMULINKového modelu dovoluje příkaz MATLABu sim.

Existují dvě varianty používání této funkce. Starší varianta 1 má syntaxi [T,X,Y] = sim('model',Timespan, Options, UT)

kde

T X Y model Timespan options UT

vektor časů výpočtu výstupních hodnot matice či struktura stavů matice či struktura výstupů název souboru obsahující model (přípona .mdl) rozsah časů simulace struktura s parametry simulace případný vstupní signál (musí obsahovat i informaci o čase)

Od verze SIMULINKu 7.3 je změněna syntaxe volání příkazu plně využívající možností objektového přístupu. Blíže viz help sim. Tento způsob dovoluje paralelní spuštění výpočtu více modelů. Nový základní tvar je SimOut = sim('model', Parameters)

kde

SimOut model Parameters

objekt obsahující výstupní informace název souboru obsahující model (přípona .mdl) a) seznam dvojic Název parametru a Hodnota b) struktura obsahující nastavení parametrů c) konfigurační sada

Využití volání SIMULINKového modelu ukažme na příkladu. Mějme v souboru zaznamenaný časový průběh vstupního a výstupního signálu (s aditivním šu‐ mem) nějakého dynamického systému. V souboru jsou na každém řádku tři hodnoty čas tk, hodnota vstupu u(tk) a výstupu y(tk). Víme, že chování systému je popsáno známým matematickým modelem

y′′′ + a2 y′′ + a1 y′ + a0 y = b2u′′ + b1u′ + b0u y′′(0) = y2 y′(0) = y1 y (0) = y0 u′(0)= u1 u (0) = u0 Máme v jednom grafu porovnat průběh zaznamenaného (vektor ym) a z modelu vypočítaného průběhu y na základě známého průběhu vstupů (vektor u). Pro výpočet výstupu systému y využijeme volání SIMULINKového modelu. SIMULINKový model Využijeme model SIMULINKu (viz Obrázek 13‐1), který představuje zápis zadané diferenciální rovnice se symbolickými názvy parametrů (včetně počátečních podmínek). Model musí obsahovat označení vstupních a výstupních signálů. Protože model obsahuje derivace vstupního signálu je vhodné pro zápis do

1

Která je podporována z důvodů zpětné kompatibility

KŘP/IMSW Modelování ve výpočtových software

13‐2 (17) 17.12.11


SIMULINKového modelu provést postupné integrace včetně přepočtu počáteč‐ ních podmínek příslušných dílčích integrátorů

y′′′ = b2u′′ − a2 y′′ + b1u′ − a1 y′ + b0u − a0 y t

y′′ = b2u′ − a2 y′ + b1u − a1 y + ∫ (b0u − a0 y )dt = b2u′ − a2 y′ + b1u − a1 y + r (t )

0

y2 = b2u1 − a2 y1 + b1u0 − a1 y0 + r (0) ⇒ r (0) = y2 + a2 y1 + a1 y0 − b2u1 − b1u0 t

y′ = b2u − a2 y + ∫ (b1u − a1 y + r )dt = b2u − a2 y + s (t ) 0

y1 = b2u0 − a2 y0 + s (t ) ⇒ s (0) = y1 + a2 y0 − b2u0 t

y = ∫ (b2u − a2 y + s (t ) )dt = 0

t t ⎛ ⎛ ⎞⎞ ⎜ ⎜ = ∫ b2u − a2 y + ∫ ⎜ b1u − a1 y + ∫ (b0u − a0 y )dt ⎟⎟ ⎟dt ⎜ ⎟ 0⎝ 0⎝ 0 ⎠⎠ t

1 u(t) b0

b1

Init. cond,: y +a y +a y -b u -b u 2

2 1

1 0

1 s

2 1

b2

Init. cond,: y +a y -b u

1 0

1

r(t)

Integrator2 a0

2 0

2 0

1 s

s(t)

Init. cond,: y 0

1 s

Integrator1 a1

y (t)

1 y(t)

Integrator a2

Obrázek 13‐1 Model pro volání z MATLABu

Starší varianta volání funkce sim Pro volání modelu mimo prostředí SIMULINKu je možné změnit některé infor‐ mace standardně zadávané při tvorbě modelu v položce hlavního menu Simulation->Configuration Parameters …. Pro tyto účely je určena speciál‐ ní struktura. Proměnnou s touto strukturou naplněnou původními hodnotami modelu zadanému jako parametr vrací funkce simget se syntaxí Options=simget(´model´)

Názvy položek struktury lze zjistit po jejím vytvoření pomocí funkce simget. Položky pak lze měnit pomocí funkce simset nebo přímým zápisem do polož‐ ky. V následujícím skriptu je ukázáno řešení zadaného příkladu využívající starší syntaxi. Je využito převzetí symbolických parametrů bloků a vstupního signálu z pracovního prostoru MATLABu. Výstupní signály jsou uvedeny jako výstupní parametry funkce sim. Doba simulace je jeden ze vstupních parametrů funkce. Kromě parametrů modelu je potřeba zadat i počáteční podmínky, které se


13‐3 (17) 17.12.11


odhadnou z dat. Výsledný graf (společný pro oba způsoby volání funkce sim) je na Obrázku 13‐2.

% skript pro příklad volání SIMULINKového modelu z MATLABu % starší varianta volání funkce sim %% vytvoření parametrů mdl='Prikl_13_1'; % název souboru s modelem (.mdl) conf=simget(mdl); % vstupní i výstupní proměnné jsou v workspace MATLABu conf=simset(conf,'SrcWorkspace','base'); conf.DstWorkspace='base'); % přímý zápis do položky conf.MaxStep=0.1; %% příprava vstupních parametrů load data.txt % načtení souboru s daty cas=data(:,1); % sloupcový vektor časů u=data(:,2); % sloupcový vektor hodnot vstupů ym=data(:,3); % sloupcový vektor hodnot zad. výstupů UU=[cas,u]; % parametry modelu a0=1.4755; a1=4.5596; a2=0.7333; b0=1.4755; b1=0.0000; b2=-5.9020; % odhady poč. podmínek z dat y0=ym(1); % odhad y(0) y1=(ym(2)-ym(1))/(cas(2)-cas(1)); % odhad první derivace y´(0) y2=0; % odhad y´´(0) u0=u(1); % odhad u(0) u1=(u(2)-u(1))/(cas(2)-cas(1)); % odhad u´(0) % starší způsob volání SIMULINKového modelu s nastavenými parametry [t,X,y] = sim(mdl,[0,60], conf, UU); % zpracování výstupních dat plot(cas,ym,t,y,'LineWidth',2)

Nová varianta volání funkce sim I v nové variantě volání funkce sim je potřeba připravit doplňující informace. Aby bylo možné model spouštět mimo prostředí SIMULINKu a s jinými parame‐ try než byly použity při jeho tvorbě, je nutné opět vytvořit konfigurační sadu s hodnotami konfiguračních parametrů. Zatímco ve staré variantě byly para‐ metry uloženy ve struktuře neobsahující všechny parametry konfigurace mode‐ lu, v nové variantě jde o objekt osahující parametry všechny. Parametry konfi‐ gurace modelu lze nastavit příkazy MATLABu podobně, jako se nastavují v položce hlavního menu Simulation->Configuration Parameters … okna modelu. Konkrétní názvy parametrů, které je potřeba použít v příkazu pro nastavení parametrů modelu se syntaxí

set_param(model,‘Nazev Parametru‘,‘Hodnota‘)

se zjistí např. kliknutím pravým tlačítkem myši nad příslušnou položkou dialo‐ gového okna Configuration Parameters … a potvrzením dotazu What‘s This?. Je potřeba nastavit aby se jako vstupní signál použily připravené pro‐ měnné v pracovním prostoru a aby se výsledek uložil do zadané výstupní pro‐ měnné (typu objekt) v požadovaných časech. Tyto a další příkazy pro vytvoření požadované konfigurace pro spuštění modelu jsou okomentované uvedeny ve výpisu skriptu pro řešení příkladu. Část zabývající se přípravou vstupních dat je stejná jako v předchozím případě.


13‐4 (17) 17.12.11


% skript pro příklad volání SIMULINKového modelu z MATLABu % nová varianta volání funkce sim %% vytvoření konfigurační sady mdl='Prikl_13_1'; % název souboru s modelem (.mdl) load_system(mdl) % načti model do paměti % uložení aktuální načtené konfigurace pro následující změny conf=getActiveConfigSet(mdl); %% nastavení nové konfigurační sady "conf" % povolení čtení vstupního signálu (z workspace) set_param(conf,'LoadExternalInput','on') % určení proměnné pro vstup (včetně času) set_param(conf,'ExternalInput','[cas,u]') % povolení zápisu výstupu (do workspace) set_param(conf,'SaveOutput','on') % určení proměnné pro uložení výstupu set_param(conf,'OutputSaveName','VYP') % určení formátu uložení výstupu set_param(conf,'SaveFormat','StructureWithTime') % uložení výstupu pouze v požadovaných časech set_param(conf,'OutputOption','SpecifiedOutputTimes') % určení požadovaných časů set_param(conf,'OutputTimes','cas') % nastavení koncového času výpočtu set_param(conf,'StopTime','cas(end)') % nastavení max. kroku výpočtu set_param(conf,'MaxStep','0.1') % nastavení Single-Output Modu set_param(conf,'ReturnWorkspaceOutputs','on') %% příprava vstupních parametrů load data.txt % načtení souboru s daty cas=data(:,1); % sloupcový vektor časů u=data(:,2); % sloupcový vektor hodnot vstupů ym=data(:,3); % sloupcový vektor hodnot měř. výstupů % parametry modelu a0=1.4755; a1=4.5596; a2=0.7333; b0=1.4755; b1=0.0000; b2=-5.9020; % odhady poč. podmínek z dat y0=ym(1); % odhad y(0) y1=(ym(2)-ym(1))/(cas(2)-cas(1)); % odhad první derivace y´(0) y2=0; % odhad y´´(0) u0=u(1); % odhad u(0) u1=(u(2)-u(1))/(cas(2)-cas(1)); % odhad u´(0) %% volání SIMULINKového modelu s nastavenými parametry Out=sim(mdl,conf); % zpracování výstupního objektu vyp=get(Out,'VYP'); % vyjmi výsledky v zadané vlastnosti y=vyp.signals.values(:,1); % výsledky ve standardní struktuře t=vyp.time; % čas výsledků (shodný s cas) plot(cas,ym,t,y,'LineWidth',2)

Výsledný graf (stejný pro obě varianty) je na Obrázku 13‐2. Je vidět, že největší neshoda je na začátku průběhů daná zejména nepřesnými hodnotami počá‐ tečních podmínek výpočtu diferenciální rovnice získanými odhadem z vstup‐ ních dat.


13‐5 (17) 17.12.11


15

10

5

0

-5

-10

-15

0

10

20

30

40

50

60

Obrázek 13‐2 Volání modelu SIMULINKu z MATLABu

13.2 Použití Symbolic Math Toolbox Symbolické výpočty nejsou zahrnuty do základního MATLABu a jsou podporo‐ vány pomocí Symbolic Math Toolbox obsahujícího samostatné výpočetní jádro. Původně (do verze v4.9 tj. R2007b) Symbolic Math Tbx obsahoval výpočetní jádro MAPLE fy Waterloo Maple, Inc., nyní (od verze 5.0 R2008a) využívá výpo‐ četní jádro MuPAD (Multiprocessing Algrebraic Data Tool) fy SciFace Software GmbH & Co. Kromě možnosti používat volání funkcí podle konvencí MATLABu MuPAD obsahuje i vlastní jazyk včetně editoru (MuPAD Notebook Interface) a možnosti konvertovat výsledky do formátu MathML a TeX pro potřeby publi‐ help symbolic či kování. Blíže viz http://www.mathworks.com/access/helpdesk/help/pdf_doc/symbo lic/mupad_tutorial.pdf

Symbolic Math Tbx podporuje výpočty zejména v následujících oblastech: ‐ numerické výpočty s proměnnou přesností (VPA variable precision ari‐ thmetic) ‐ symbolické operace lineární algebry ‐ řešení soustav algebraických rovnic ‐ integrální transformace (Laplaceova, Fourierova, Z‐transformace atd.) ‐ řešení integrálů ‐ řešení diferenciálních rovnic ‐ zjednodušení symbolických výrazů Co znamenají výpočty s požadovanou přesností? Standardní výpočty se v MATLABu provádějí v tzv. double precision využívající pro uložení 64 bitů. Při těchto výpočtech je přesnost cca 15 platných cifer. V případě výpočtu s požadovanou přesností jde o numerické výpočty prováděné na zadaný počet


13‐6 (17) 17.12.11


platných cifer. Vyhodnocení výrazu s požadovanou přesností je možné provést pomocí funkce vpa se syntaxí

v=vpa(S,D)

kde

S je vyhodnocovaný výraz D je nepovinný počet požadovaných cifer; je‐li vynechán, použije se na‐ stavená hodnota. Hodnota se nastavuje či zjišťuje funkcí digits. Blíže viz help vpa.

Je potřeba dát si pozor na to, aby výpočet výrazu byl prováděn s proměnnými typu sym a ne standardního typu double. V druhém případě je převeden do požadované přesnosti jen výsledek výpočtu což není přesný výsledek. Ukažme

(

)

tento problém na příkladě výpočtu hodnoty výrazu 1+ 5 2 . Nastavme vý‐ pis na maximální počet číslic a proveďme standardní výpočet MATLABu (datový typ double přesnost cca 15 číslic) >> format long >> v=(1+sqrt(5))/2 v = 1.618033988749895

standardní výpočet v double s převodem výsledku do implicitní přesnosti vý‐ počtů s proměnnou přesností (32) = nepřesný výsledek od 17‐té číslice >> v1=vpa((1+sqrt(5))/2) v1 = 1.6180339887498949025257388711907

správný výpočet s implicitní přesností 32 číslic – interpretace zápisu (řetězec) v MuPADu >> v2=vpa('(1+sqrt(5))/2') v2 = 1.6180339887498948482045868343656

jiný správný výpočet s požadovanou přesností 48 číslic – definice čísla 5 jako datového typu sym a výpočet obsahující proměnnou typu sym se provede v MuPADu >> s=sym(5); >> v3=vpa((1+sqrt(s))/2,48) v3 = 1.61803398874989484820458683436563811772030917981

Co se rozumí pod pojmem symbolický výpočet? Jednoduše řečeno při symbo‐ lickém výpočtu výsledkem výrazu (vzorce) není číslo ale jiný výraz (vzorec). Např. chceme‐li určit výsledek výrazu a^2+2*a*b+b^2 tak v případě numeric‐ kého výpočtu musíme mít proměnné a a b naplněny čísly a číslo také jako vý‐ sledek dostaneme. V případě symbolického výpočtu dostaneme jako výsledek výraz (a+b)^2. Aby bylo možné používat standardní matematické operace a názvy funkcí pro provádění symbolických výpočtů, musí být proměnné typu symbolický objekt. Vytvoření symbolických proměnných Jsou dvě možnosti definovat proměnnou jako symbolický objekt, buď jednotli‐ vě pomocí funkce sym nebo hromadně pomocí příkazu syms. Blíže viz help sym a help syms >> x=sym('x'); y= sym('y'); z=sym('z'); >> syms a b c


13‐7 (17) 17.12.11


Zjednodušování a výpis symbolických výrazů Obvykle je snaha používat co nejjednodušší vyjádření výrazů. Pro zjednodušení je možné aplikovat několik postupů. Tbx nabízí několik funkcí pro zjednodušení (simplify, factor, combine, expand atd.) ale pokud nevíme, kterou použít, je nejjednodušší použít funkci simple, které vyzkouší všechny možnos‐ ti. Určitým nedostatkem 2 MATLABu vyplývajícím z textové formy komunikace je ne příliš přehledná textová forma zadávání i výpisu výrazů. Funkce pretty zobrazuje (byť opět textově) výraz poněkud přehledněji. Blíže viz help simple a help pretty. Použití obou funkcí ukažme na příkladě zobrazení a

(a + b )2 ⎛ zjednodušení výrazu

−1

a+b ⎞ ⎟ 2 2 ⎜ a − b ⎝ ab − b ⎠

>> v=sym('(a+b)^2/(a^2-b^2)*((a+b)/(a*b-b^2))^(-1)') v = ((a + b)*(- b^2 + a*b))/(a^2 - b^2) >> pretty(v) 2 (a + b) (- b + a b) -------------------2 2 a - b >> vysl=simple(v) simplify: b radsimp: ((a + b)*(- b^2 + a*b))/(a^2 - b^2) simplify(100): b combine(sincos): ((a + b)*(- b^2 + a*b))/(a^2 - b^2) combine(sinhcosh): ((a + b)*(- b^2 + a*b))/(a^2 - b^2) combine(ln): ((a + b)*(- b^2 + a*b))/(a^2 - b^2) factor: b expand: (a^2*b)/(a^2 - b^2) - b^3/(a^2 - b^2) combine: ((a + b)*(- b^2 + a*b))/(a^2 - b^2) rewrite(exp): ((a + b)*(- b^2 + a*b))/(a^2 - b^2) rewrite(sincos): ((a + b)*(- b^2 + a*b))/(a^2 - b^2) rewrite(sinhcosh): ((a + b)*(- b^2 + a*b))/(a^2 - b^2) rewrite(tan): ((a + b)*(- b^2 + a*b))/(a^2 - b^2) mwcos2sin: ((a + b)*(- b^2 + a*b))/(a^2 - b^2) collect(b): b vysl = b

Je možné také vyjádřit symbolický výraz ve formě syntaxe jazyka „C“ (funkce ccode) nebo převést do vyjádření pro systém LaTex (funkce latex).

13.2.1 Lineární algebra a algebraické rovnice Možnosti symbolických výpočtů v oblasti lineární algebry a řešení rovnic si uká‐ žeme na několika příkladech. V těchto případech nejsou potřeba zvláštní funk‐ ce – pouze je potřeba, aby se standardní operace či funkce jako je např. de‐ terminant (funkce det), inverze (funkce inv), vlastní čísla (funkce eig) a cha‐ rakteristický polynom (funkce poly) prováděly nad proměnnými typu sym. Je možné definovat matici symbolických proměnných a dále s ní pracovat. 2

Tato nevýhoda je výhodou v případě, že potřebujeme výraz použít pro výpočet v programu


13‐8 (17) 17.12.11


Standardní výpočty se symbolickými proměnnými vytvoření matice a vektoru symbolických proměnných >> syms a11 a12 a21 a22 b1 b2 >> A=[a11,a12;a21,a22]; b=[b1;b2];

násobení matic >> c=A*b c = a11*b1 + a21*b1 + >> c=A*[2;5] c = 2*a11 + 2*a21 +

a12*b2 a22*b2 5*a12 5*a22

výpočet determinantu >> M=sym('[a b;x y]') M = [ a, b] [ x, y] >> D=det(M) D = a*y - b*x

výpočet inverzní matice >> I=inv(M) I = [ y/(a*y - b*x), -b/(a*y - b*x)] [ -x/(a*y - b*x), a/(a*y - b*x)]

výpočet vlastních čísel matice >> E=eig(M) E = a/2 + y/2 - (a^2 - 2*a*y + y^2 + 4*b*x)^(1/2)/2 a/2 + y/2 + (a^2 - 2*a*y + y^2 + 4*b*x)^(1/2)/2

výpočet charakteristického polynomu matice >> P=poly(M) P = t^2 + (- a - y)*t + a*y - b*x

Algebraické rovnice Důležitou oblastí je řešení algebraických rovnic. Pokud řešení existuje, lze je získat ve formě symbolického výrazu (obvykle velmi složitého) pomocí funkce solve se syntaxí r=solve(eq1,eq2,…,eqn,var1,var2,…,varn) kde eq jsou definice rovnic a var označení proměnných, blíže viz help solve. Použití ukážeme opět na několika příkladech. Úloha je‐li známa jedna rovnice f(p,x)=0 či soustava rovnic f(p,x)=0 nalezněte všechna analytická řešení x=g(p) Příklad nalezněte kořeny obecné kvadratické rovnice ax 2 + bx + c = 0 (dvě řešení) Řešení definice symbolických proměnných >> syms a b c

výpočet kvadratické rovnice ax 2 + bx + c = 0 >> x2=solve('a*x^2+b*x+c=0','x') x2 = -(b + (b^2 - 4*a*c)^(1/2))/(2*a) -(b - (b^2 - 4*a*c)^(1/2))/(2*a)


13‐9 (17) 17.12.11


Příklad nalezněte kořeny kubické rovnice x + ax + bx + c = 0 (tři řešení) 3

2

Řešení definice symbolických proměnných >> syms a b c

výpočet kubické rovnice x 3 + ax 2 + bx + c = 0 >> x3=solve('x^3+a*x^2+b*x+c=0','x') x3 = ((a*b)/6 - c/2 + ((b/3 - a^2/9)^3 + (a^3/27 - (b*a)/6 + c/2)^2)^(1/2) - a^3/27)^(1/3) - (b/3 - a^2/9)/((a*b)/6 - c/2 + ((b/3 - a^2/9)^3 + (a^3/27 - (b*a)/6 + c/2)^2)^(1/2) a^3/27)^(1/3) - a/3 (b/3 - a^2/9)/(2*((a*b)/6 - c/2 + ((b/3 - a^2/9)^3 + (a^3/27 - (b*a)/6 + c/2)^2)^(1/2) - a^3/27)^(1/3)) - a/3 ((a*b)/6 - c/2 + ((b/3 - a^2/9)^3 + (a^3/27 - (b*a)/6 + c/2)^2)^(1/2) a^3/27)^(1/3)/2 (3^(1/2)*((b/3 a^2/9)/((a*b)/6 - c/2 + ((b/3 - a^2/9)^3 + (a^3/27 - (b*a)/6 + c/2)^2)^(1/2) - a^3/27)^(1/3) + ((a*b)/6 - c/2 + ((b/3 - a^2/9)^3 + (a^3/27 - (b*a)/6 + c/2)^2)^(1/2) - a^3/27)^(1/3))*i)/2 (b/3 - a^2/9)/(2*((a*b)/6 - c/2 + ((b/3 - a^2/9)^3 + (a^3/27 - (b*a)/6 + c/2)^2)^(1/2) - a^3/27)^(1/3)) - a/3 ((a*b)/6 - c/2 + ((b/3 - a^2/9)^3 + (a^3/27 - (b*a)/6 + c/2)^2)^(1/2) a^3/27)^(1/3)/2 + (3^(1/2)*((b/3 a^2/9)/((a*b)/6 - c/2 + ((b/3 - a^2/9)^3 + (a^3/27 - (b*a)/6 + c/2)^2)^(1/2) - a^3/27)^(1/3) + ((a*b)/6 - c/2 + ((b/3 - a^2/9)^3 + (a^3/27 - (b*a)/6 + c/2)^2)^(1/2) - a^3/27)^(1/3))*i)/2

Příklad nalezněte kořeny kubické rovnice x 3 + 7 x 2 − x / 5 + 11 = 0 (tři přesná numerická řešení ) Řešení >> v=solve('x^3+7*x^2-x/5+11=0'); v = - 248/(45*(2489/135 (3375^(1/2)*582319^(1/2))/3375)^(1/3)) - (2489/135 (3375^(1/2)*582319^(1/2))/3375)^(1/3) - 7/3 (3^(1/2)*(248/(45*(2489/135 (3375^(1/2)*582319^(1/2))/3375)^(1/3)) - (2489/135 (3375^(1/2)*582319^(1/2))/3375)^(1/3))*i)/2 + 124/(45*(2489/135 (3375^(1/2)*582319^(1/2))/3375)^(1/3)) + (2489/135 (3375^(1/2)*582319^(1/2))/3375)^(1/3)/2 - 7/3 124/(45*(2489/135 (3375^(1/2)*582319^(1/2))/3375)^(1/3)) (3^(1/2)*(248/(45*(2489/135 (3375^(1/2)*582319^(1/2))/3375)^(1/3)) - (2489/135 (3375^(1/2)*582319^(1/2))/3375)^(1/3))*i)/2 + (2489/135 (3375^(1/2)*582319^(1/2))/3375)^(1/3)/2 - 7/3

což převedeno do standardního číselného vyjádření funkcí double je přibližně >> double(v) ans = -7.2376 0.1188 + 1.2271i 0.1188 - 1.2271i

Příklad nalezněte všechna řešení soustavy dvou nelineárních rovnic x 2 + ay 2 = 1 y 2 − bx 2 = 1 (čtyři páry řešení) Řešení definice symbolických proměnných >> syms a b 2 2 výpočet soustavy rovnic x + ay = 1

y 2 − bx 2 = 1

>> eq1=sym('x^2+a*y^2=1'); >> eq2=sym('y^2-b*x^2=1'); >> v=solve(eq1,eq2,'x','y') v = x: [4x1 sym]


13‐10 (17) 17.12.11


y: [4x1 sym] >> v.x ans = (-(a - 1)/(a*b -(-(a - 1)/(a*b (-(a - 1)/(a*b -(-(a - 1)/(a*b

+ + + +

>> v.y 1))^(1/2) ans = ((b + 1)/(a*b 1))^(1/2) ((b + 1)/(a*b 1))^(1/2) -((b + 1)/(a*b 1))^(1/2) -((b + 1)/(a*b

Příklad nalezněte všechna řešení soustavy dvou rovnic

+ + + +

1))^(1/2) 1))^(1/2) 1))^(1/2) 1))^(1/2)

(1 + ax )y = c (1 + by )x = c

Řešení >> v=solve ('(1+a*x)*y=c','(1+b*y)*x=c') v = x: [2x1 sym] y: [2x1 sym] >> v.x ans = -((b*c)/2 - (a*c)/2 + (a^2*c^2 - 2*a*b*c^2 + 2*a*c + b^2*c^2 + 2*b*c + 1)^(1/2)/2 + 1/2)/a ((a*c)/2 - (b*c)/2 + (a^2*c^2 - 2*a*b*c^2 + 2*a*c + b^2*c^2 + 2*b*c + 1)^(1/2)/2 - 1/2)/a >> v.y ans = -(a*c - b*c + (a^2*c^2 - 2*a*b*c^2 + 2*a*c + b^2*c^2 + 2*b*c + 1)^(1/2) + 1)/(2*b) -(a*c - b*c - (a^2*c^2 - 2*a*b*c^2 + 2*a*c + b^2*c^2 + 2*b*c + 1)^(1/2) + 1)/(2*b)

13.2.2 Derivace, integrály, integrální transformace Pomocí Symbolic Math Tbx lze také určit symbolické (analytické) derivace a integrace funkcí. Blíže viz help diff a help int. Určení derivace funkce se provádí pomocí funkce diff r=diff(S,var,n) kde S je symbolický výraz, který se má derivovat, volitelný parametr var je označení proměnné podle které se má derivovat a volitelný parametr 1 ≤ n je řád derivace. Úloha je‐li známa funkce f(p,x), nalezněte n‐tou derivaci podle zvolené proměnné Příklad mějme funkci f (a, x ) = cos (ax ) , určete první derivaci podle x, dru‐ hou derivaci podle x a první derivaci podle a a

Řešení první derivace podle x >> dx1=diff(sym('cos(a*x)â')) dx1 = -a^2*cos(a*x)^(a - 1)*sin(a*x)

df dx = −a 2 cosa−1 (ax )sin (ax )

druhá derivace podle x >> dx2=diff(sym('cos(a*x)â'),2) dx2 = a^3*cos(a*x)^(a 2)*sin(a*x)^2*(a a^3*cos(a*x)*cos(a*x)^(a - 1)

-

1)

-

d 2 f dx2 = a3 cosa−2 (ax)sin(ax)(a −1) − a3 cosa−2 (ax) cosa−1 (ax)

první derivace podle a >> da1=diff(sym('cos(a*x)â'),'a') da1 = log(cos(a*x))*cos(a*x)â - a*x*cos(a*x)^(a - 1)*sin(a*x) >> pretty(da1)


13‐11 (17) 17.12.11


a a - 1 log(cos(a x)) cos(a x) - a x cos(a x) sin(a x)

(

)

df da = ln cosa−1 (ax ) cosa − ax cosa−1 (ax )sin (ax )

Určení integrálu funkce se provádí pomocí funkce int se syntaxí r=int(S,var,a,b) kde S je symbolický výraz, který se má integrovat, volitelný parametr var je označení proměnné, podle které se má integrovat a volitelné parametry a, b jsou dolní a horní integrační mez v případě, že chceme určitý integrál. Úloha je‐li známa funkce f(p,x), nalezněte neurčitý či určitý integrál podle proměnné x

(

∫ sin (a + 2

)

f (a, x ) = sin 2 a + x , určete neurčitý integrál

Příklad mějme funkci

)

x dx

Řešení definice definice symbolických proměnných >> syms a x >> in=int(sym('sin(a+sqrt(x))^2'),'x') in = x/2 - (x^(1/2)*sin(2*a + 2*x^(1/2)))/2 + sin(a + x^(1/2))^2/2 >> pretty(in) 1/2 1/2 1/2 2 x x sin(2 a + 2 x ) sin(a + x ) - - ---------------------- + -------------2 2 2

∫ sin (a + 2

)

x dx =

[

(

)

Příklad mějme funkci f (R, x ) =

)]

R 2 − x 2 , určete hodnotu určitého integrálu

R

∫

(

1 x − x sin 2a + 2 x + sin 2 a + x 2

x

R 2 − x 2 dx a integrálu jako funkce horní meze

−R

∫

R 2 − t 2 dt

−R

Řešení definice symbolických proměnných >> syms R x >> eq=sym('sqrt(R^2-x^2)'); >> in1=int(eq,'x',-R,R) in1 = (pi*R^2)/2 R

∫

−R

1 R 2 − x 2 dx = πR 3 2

>> in2=int(eq,'x',-R,x) in2 = (R^2*asin(R/(R^2)^(1/2)))/2 + (R^2*asin(x/(R^2)^(1/2)))/2 (x*(R^2 - x^2)^(1/2))/2 >> pretty(in2) 2 / R \ 2 / x \ R asin| ------- | R asin| ------- | | 2 1/2 | | 2 1/2 | 2 2 1/2 \ (R ) / \ (R ) / x (R - x ) ------------------ + ------------------ + -------------2 2 2


13‐12 (17) 17.12.11

+


x

∫

R 2 − t 2 dt =

−R

⎤ ⎛ R ⎞ ⎛ x ⎞ 1⎡ 2 ⎟⎟ + R 2 arcsin ⎜⎜ ⎟⎟ + x R 2 − x 2 ⎥ ⎢ R arcsin ⎜⎜ 2 2 2 ⎢⎣ ⎥⎦ ⎝ R ⎠ ⎝ R ⎠

Integrální transformace Symbolic Math Tbx obsahuje funkce pro výpočet normálních i zpětných inte‐ grálních transformací – Fourierovi, Laplaceovi a Z‐transformace. Blíže viz help fourier, help ifourier, help laplace, help ilaplace, help ztrans a help iztrans Fourierova transformace je realizována funkcí F=fourier(f) a zpětná Fourierova transformace je realizována funkcí f=ifourier(F) kde f(x) je původní funkce a F(w) je Fourierův obraz Úloha je‐li známa funkce f(x), nalezněte Fourierův obraz F(w) podle vztahu

F (w ) =

+∞

∫ f (x )e

−iwx

dx a je‐li znám obraz F(w) nalezněte originál f(x) podle

−∞

vztahu f ( x ) =

1 2π

+∞

∫ F (w)e

iwx

dw

−∞

2 −x Příklad mějme funkci f ( x ) = x e , nalezněte Fourierův obraz F(w) a k tomu‐ to obrazu nalezněte originál 2

Řešení definice symbolických proměnných >> syms x >> f=x^2*exp(-x^2); >> F=fourier(f) F = pi^(1/2)/(2*exp(w^2/4)) - (pi^(1/2)*w^2)/(4*exp(w^2/4))

f (x ) = x e

2 − x2

⇒

w2

π⎛

w2 ⎞ − 4 ⎜⎜1 − ⎟e F (w ) = 2 ⎝ 2 ⎟⎠

>> fo=ifourier(F) fo = (pi/exp(x^2) (pi^(1/2)*((4*pi^(1/2))/exp(x^2) (8*pi^(1/2)*x^2)/exp(x^2)))/4)/(2*pi) >> fo=simple(fo) fo = x^2/exp(x^2)

-

Laplaceova transformace je realizována funkcí L=laplace(f) a zpětná Laplaceova transformace je realizována funkcí f=ilaplace(L) kde f(t) je původní funkce a L(s) je Laplaceův obraz


13‐13 (17) 17.12.11


Úloha je‐li známa funkce f(t), nalezněte Laplaceův obraz L(s) podle vztahu

L (s ) =

+∞

∫ f (t )e

− st

dt je‐li znám obraz L(s) nalezněte originál f(t) podle vztahu

0

c + i∞

1 f (t ) = L(s )e st ds ∫ 2πi c − i∞

2 − at Příklad mějme funkci f (t ) = t e , nalezněte Laplaceův obraz L(s) a k tomuto obrazu nalezněte originál

Řešení definice symbolických proměnných >> syms a t >> f=t^2*exp(-a*t); >> L=laplace(f) L = 2/(a + s)^3 >> fo=ilaplace(L) fo = t^2/exp(a*t)

Z‐ transformace je realizována funkcí F=ztrans(f) a zpětná Z‐transformace je realizována funkcí f=iztrans(F) kde f(n) je původní funkce a F(z) je obraz v Z‐transformaci Úloha je‐li známa funkce f(n), nalezněte obraz Z‐transformaci F(z) podle vztahu F ( z ) =

∞

∑ f (n )z

−n

a je‐li znám obraz F(z) nalezněte originál f(n)

n=0

podle vztahu f (n ) =

1 2πi

∫ F (z )z

n −1

z =R

dz n = 1,2,K

Příklad mějme funkci f (n) = n + an , nalezněte obraz v Z‐transformaci F(z) a k tomuto obrazu nalezněte originál 2

Řešení definice symbolických proměnných >> syms n a >> f=n^2+a*n; >> F=ztrans(f) F = (z^2 + z)/(z - 1)^3 + (a*z)/(z - 1)^2 >> fo=iztrans(F) fo = 2*binomial(n-1,2)-(kroneckerDelta(n,0)-1)*(a+1)+ (a+3)*(n+kroneckerDelta(n,0)-1)-2*kroneckerDelta(n,0) >> fo=simple(fo) fo = n*(a + n)

13.2.3 Diferenciální rovnice Poslední funkcí Symbolic Math Tbx, o které se zmíníme, je funkce dsolve pro řešení diferenciálních rovnic. f=dsolve(eq1,eq2,...,cond1,cond2,...,var)


13‐14 (17) 17.12.11


kde eq jsou definice rovnic, cond definice počátečních či okrajových podmínek a var označení nezávisle proměnné, blíže viz help dsolve.

Úloha je‐li známa diferenciální rovnice f(x,y,y΄,..)=0 nebo soustava diferenciálních rovnic případně s počátečními podmínkami nalezněte funkce vyhovující daným rovnicím Příklad mějme rovnici y′′ + (1 − y′) = 1 , nalezněte obecné řešení (dvě řešení) 2

Řešení >> r=dsolve('D2y+(1-Dy)^2=1') r = C12 log(C11 - (C10*exp(2*t))/2)

y1 (t ) = C

(

)

y2 (t ) = ln C1 − C2e 2t z′ = a − y , nalezněte obecné řešení

Příklad mějme rovnice y′ = az

Řešení definice symbolických proměnných >> syms a >> e1='Dy=a*z'; e2='Dz=y-a'; >> r=dsolve(e1,e2) r = z: [1x1 sym] y: [1x1 sym] >> z=r.z z = C59*exp(a^(1/2)*t) + C60/exp(a^(1/2)*t)

z (t ) = C1et

a

+ C2 e − t

a

>> y=r.y y = a + C59*a^(1/2)*exp(a^(1/2)*t) - (C60*a^(1/2))/exp(a^(1/2)*t)

y (t ) = a + C1 a et

a

+ C2 a e − t

a

Příklad mějme rovnice y′ = az z′ = a − y a počáteční podmínky y(0)=2, z(0)=‐3. Nalezněte řešení pro nezávisle proměnnou x. Řešení definice symbolických proměnných >> syms a >> e1='Dy=a*z'; e2='Dz=y-a'; >> r=dsolve(e1,e2,'y(0)=2','z(0)=-3','x') r = z: [1x1 sym] y: [1x1 sym] >> z=r.z z = -(3*a^(1/2)-a+2)/(2*a^(1/2)*exp(a^(1/2)*x)) (exp(a^(1/2)*x)*(a+3*a^(1/2)-2))/(2*a^(1/2))

z (t ) = −

3 a +a−2 x e 2 a

a

−

3 a − a + 2 −x e 2 a

-

a

>> y=r.y y = a-(exp(a^(1/2)*x)*(a+3*a^(1/2)-2))/2+(3*a^(1/2)-a+2) /(2*exp(a^(1/2)*x))

y (t ) = a −

3 a +a−2 x e 2

a

+

3 a − a + 2 −x e 2


a

13‐15 (17) 17.12.11


dy =axy y (t = 0 ) = c dt . Nalezněte řešení pro Příklad mějme rovnice dx = a ( x + y ) x (t = 0 ) = c dt nezávisle proměnné x a y. Řešení >> v=dsolve('Dy=a*(y-x)','Dx=a*(x+y)','y(0)=c','x(0)=c') v = x: [1x1 sym] y: [1x1 sym] >> v.x ans = c*exp(a*t*(1-i))*(1/2+i/2) + c*exp(a*t*(1+i))*(1/2-i/2) >> x=simple(v.x) x = c*exp(a*t)*(cos(a*t) + sin(a*t)) >> v.y ans = c*exp(a*t*(1-i))*(1/2-i/2) + c*exp(a*t*(1+i))*(1/2+i/2) >> y=simple(v.y) y = c*exp(a*t)*(cos(a*t) - sin(a*t))


13‐16 (17) 17.12.11


Pojmy k zapamatování Okruhy problémů: použití modelu SIMULINKu z MATLABu, symbolické výpočty v Symbolic Math Tbx Použité nástroje:

Funkce: diff, digits, double, dsolve, int, pretty, set_param, sim, simget, simple, simset, solve, sym, syms, vpa

Otázky na procvičení 1. 2. 3. 4. 5. 6. 7.

Jak se nazývá funkce umožňující volání modelu SIMULINKu z příkazové řádky MATLABu? Jaké jsou výhody novější varianty volání modelu SIMULINKu z MATLABu? Co znamená výpočet s požadovanou přesností a jakou se jmenuje funkce pro nastavení přesnosti? Jaká je standardní přesnost výpočtů v MATLABu? Jakého typu jsou proměnné použité pro symbolické výpočty a proč? Jak je možné, že stejná funkce (např. inv) dá jednou numerický a jindy symbolický výsledek? Co zajišťuje funkce simple?

Odkazy a další studijní prameny • •

• •

elektronická dokumentace k programu SIMULINK – v „Functions / Simu‐ lations“ funkce „sim“ elektronická dokumentace k programu MATLAB – v „Symbolic Math Toolbox / Getting Started“ část „Symbolic Objects“ a „Creating Symbolic Variables and Expressions“ elektronická dokumentace k programu MATLAB – nápověda k jednotli‐ vým funkcím elektronická dokumentace – www.mathworks.com/help/


13‐17 (17) 17.12.11


13 Další možnosti MATLABu a SIMULINKu

Recommend Documents