Modelování 2: Nelineární algebraická rovnice

Po úvodním díle o funkcích se teď zaměříme na řešení soustavy nelineárních rovnic. V této části začneme s jednou nelineární rovnicí a soustavu rovnic si necháme na příště. Také si ukážeme a vysvětlíme základní strukturu programu, podle které budeme řešit i všechny ostatní typy úloh v předmětu “Modelování”.

Jak budeme organizovat náš kód?

Začněme s tím, jak by měla vypadat struktura programu. Představme si, že chceme najít řešení nějaké konkrétní nelineární rovnice, například $3x^3+\exp{(-x)}=0$ . V předmětu zaměřeném na numerickou matematiku bychom si nejdříve vybrali jednu z mnoha metod na hledání kořene nelineární rovnice a algoritmus této metody zapsali ve formě počítačového programu, kerý úlohu vyřeší. Určitě vás napadne, že při psaní programu je lepší co nejvice oddělit algoritmickou část týkající se dané metody od části specifické pro konkrétní nelineární rovnici. Jestliže máme fungující část kódu implementující numerickou metodu, tak ji můžeme znovu použít i později pro řešení jiných rovnic. Výhodou takového rozdělení je to, že už si nemusíme pamatovat přesné detaily algoritmu, ale stačí nám základní povědomí o dané metodě.

V současnosti už pro běžné účely nebude nutné programovat algoritmus numerické metody, protože existuje řada knihoven numerických algoritmů napsaných pro nejrůznější programovací jazyky. Tyto knihovny jsou prověřené a optimalizované desetiletími používání v nejrůznějších aplikacích. Užívateli psané programy mají od nepaměti následující strukturu:

Hlavní program, tj. část programu, který něco provádí a dostane se k okazmžiku, kdy je potřeba vyřešit nějakou numerickou úlohu, například najít kořeny nelineární rovnice f(x)=0. Tato část programu obsahuje volání funkce z numerické knihovny, která úlohu vyřeší a vrátí výsledek. Hlavní program potom může pokračovat dále a pracovat s tímto výsledkem.
Vlastní funkce z numerické knihovny (numerická funkce) implementující numerickou metodu, například Newtonovu metodu. Tuto funkci uživatel neprogramuje, pouze si nastuduje její dokumentaci, aby věděl, jak funkci používat.
“Modelová” funkce, tj. uživatelem napsaná funkce, která definuje daný numerický problém. V našem příkladě bude obsahovat konkrétní nelineární rovnici f(x)=0. Po zadání hodnoty x (vstupní proměnná) vrátí residuál f(x) (výstupní proměnná). Kdykoliv algoritmus uvnitř numerické funkce potřebuje znát hodnotu funkce v konkrétním bodě x, f(x), zavolá tuto uživatelem dodanou funkci.

Tuto strukturu budeme používat i pro naše programy v MATLABu v předmětu “Modelování”. Naučíme se, jak volat numerické funkce řešící různé typy numerických problémů a také, jak vytvořit “modelové” funkce pro konkrétní úlohy.

Jednoduchá nelineární rovnice

Ukažme si, jak najít řešení nelineární rovnice z úvodu článku. Nejdříve je nutné zapsat rovnici ve tvaru f(x)=0, pro náš příklad je to $f(x)=3x^3+\exp{(-x)}=0$

K hledání řešení použijeme numerickou funkci fzero, která je k dispozici v základní sestavě MATLABu. Vstupní parametry, které musíme funkci dodat jsou

odkaz na funkci residuálu, tj. “modelovou” funkci f(x)=0 a
počáteční odhad.

Psaní funkcí jsme vysvětlili v minulém díle. Funkci residuálu můžeme zapsat do souboru se jménem fun.m v aktuálním adresáři

function res = fun(x)
% Example of residual function for a nonlinear equation
%
res = 3*x.^3 + exp(-x);

Vidíme, že funkce musí mít jednu vstupní a jednu výstupní proměnnou. Zvolíme-li počáteční odhad x=0, hledání řešení můžeme spustit příkazem

>> fzero(@fun, 0)

Znak zavináče @ znamená, že příslušná proměnná je typu odkaz na funkci (function handle) místo normálního typu. Všimněte si, že jsme v definici funkce fun použili operátor tečkového umocňování .^ místo maticového ^. Proto můžeme také zobrazit graf naší funkce

>> x = linspace(-7,3,100); >> y = fun(x); >> plot(x,y)

Vidíme, že rovnice má více řešení. Další řešení dostaneme, když jako druhý parametr zadáme interval, ve kterém se řešení nachází.

>> fzero(@fun, [-10,-5])

Anonymní funkce

Pokud máme jednoduchou funkci, která se skládá jen z jednoho řádku a nechceme pro ni vytvářet zvláštní soubor na disku, můžeme použít koncept tzv. anonymních funkcí. Anonymní funkci nadefinujeme přímo na příkazovém řádku a odkaz na tuto funkci uložíme do proměnné fun2. Potom zavoláme funkci fzero.

>> fun2 = @(x) 3*x.^3 + exp(-x) >> fzero(fun2, 0)

Znakem zavináče @ se můžeme odkazovat na existující funkci (např. @fun). Pokud ale za znakem @ není jméno funkce, ale seznam proměnných uzavřený v závorkách, tj. @(x), potom odkazujeme na funkci, kterou budeme právě definovat. Tato funkce bude mít jednu vstupní proměnnou. Je potřeba si uvědomit, že fun2 není funkce, ale proměnná odkazující na anonymní funkci uloženou v paměti. Ve výrazech a vzorcích můžeme fun2 používat stejně, jako by se jednalo o funkci, například y=fun2(x). Všimněte si ale, že výraz ve volání fzero už musíme psát bez zavináče @ před fun2.

Výpočet bodu varu vícesložkové směsi

Budeme pokračovat s příkladem chemicko-inženýrského problému, kterým bude určení bodu varu vícesložkové směsi.

Zadání

Určete bod varu a odpovídající složení parní fáze směsi hexan/heptan/oktan. Předpokládejte atmosférický tlak a ekvimolární složení kapalné fáze.

Matematický model

Budeme předpokládat Raoultův zákon pro rovnováhu kapalina-pára

$p y_i = p_i^\circ(t) x_i\qquad i \in {1,2,3}$ .

Sečteme-li tyto tři rovnice a uplatníme-li podmínku, že součet molárních zlomků v parní fázi musí být roven jedné, dostaneme nelineární rovnici

$p = \displaystyle\sum_{i=1}^3 p_i^\circ(t) x_i$ ,

kde jedinou neznámou zůstane teplota varu t. Závislost tenze nasycených par na teplotě určíme z Antoineho rovnice

$\ln p_i^\circ(T) = A_i - {B_i}/(T+C_i)$ .

Funkce pro výpočet residuálu

Pro výpočet tenze par použijeme funkce ptension a ptension_coef, které jsme vytvořili v minulé části. Nakopírujte proto soubory

ptension.m
ptension_coef.m

do pracovního adresáře. Ve stejném adresáři teď napíšeme funkci pro výpočet residuálu res_boiling.m

function res = res_boiling(t)
%RES_BOILING  Residual function for boiling temperature calculation problem
%  Boiling temperature calculation by solving non-linear equation
%    t is temperature in Celsius
%    Raoult law for liquid/vapour equilibrium assumed
 
patm = 101.325;    % atmospheric pressure in kPa
x1 = 0.333;
x2 = 0.333;
x3 = 1 - x1 - x2;  % liquid phase mole fractions
 
[A,B,C] = ptension_coef(); % get Antoine equation coefficients 
p1 = ptension(A(1), B(1), C(1), t); % hexane vapour pressure in kPa
p2 = ptension(A(2), B(2), C(2), t); % heptane vapour pressure in kPa
p3 = ptension(A(3), B(3), C(3), t); % octane vapour pressure in kPa
 
% calculate the residual
res = p1*x1 + p2*x2 + p3*x3 - patm;
end

V této ukázce jsme kvůli zjednodušení záměrně omezili používání maticových operací a pracujeme pouze se skaláry. Až získáte více zkušeností s MATLABEM, můžete funkci residuálu napsat i takto

function res = res_boiling(t)
%RES_BOILING  Residual function for boiling temperature calculation problem
%  Boiling temperature calculation by solving non-linear equation
%    t is temperature in Celsius
%    Raoult law for liquid/vapour equilibrium assumed
% [Version using matrix operators. MATLAB beginners might find the code
%  more difficult to understand than in the previous version]
 
patm = 101.325;          % atmospheric pressure in kPa
x(1) = 0.333; 
x(2) = 0.333;
x(3) = 1 - x(1) - x(2);  % liquid phase mole fractions vector
 
[A,B,C] = ptension_coef(); % get Antoine equation coefficients 
p = ptension(A, B, C, t);  % vector of vapour pressures in kPa
res = p*x' - patm;         % calculate the residual
end

Je na vás, které verzi budete dávát přednost. První verze bude pro účely předmětu “Modelování” naprosto dostačující a je v souladu s filosofií “Keep it simple…”. Používání maticových operací ve druhé verzi ušetří pár řádek kódu, ale zhoršuje jeho čitelnost pokud s MATLABEM začínáte. Každá řádka, kterou ve svých programech napíšete, musí být pro vás naprosto srozumitelná. Proto když si nejste naprosto jistí, je lepší se komplikovanějším konstrukcím vyhýbat. Nijak vás to v tuto chvíli neomezí, předmět “Modelování” není zaměřen na výuku programování.

Hlavní program

Než zavoláme funkci fzero, je vhodné naši funkci residuálu vyzkoušet. Má jeden vstupní parametr, kterým je teplota. Pokud bychom zadali teplotu varu směsi, tj. znali řešení, pak rovnice bude splněna a residuál bude nula. Pokud při jedné teplotě t₁ bude residuál kladný a při jiné teplotě t₂ záporný, potom řešení rovnice bude ležet někde mezi oběmi teplotami. Vyzkoušejte vámi vytvořenou funkci pro několik teplot a podívejte se, zda dokážete najít interval, ve kterém má rovnice řešení

>> res_boiling(50) >> res_boiling(100)

Pokud jsme funkci residuálu otestovali a funguje bez chyb pro různé teploty, můžeme najít řešení s funkcí fzero. Ze znalosti fyzikálního chování systému určíme realistický počáteční odhad, například 50°C

>> tb = fzero(@res_boiling, 50)

Řešení jsme uložili do proměnné tb. Po použití numerické metody je potřeba zkontrolovat, zda příslušná metoda zkonvertovala ke správnému řešení a ověřit, že residuál je blízko nuly.

>> res_boiling(tb)

Zbývá určit složení parní fáze, což už je, jakmile známe teplotu varu, jednoduché. Znovu můžeme použít naše funkce pro výpočet tenze par, tentokrát voláním z hlavního programu (příkazového řádku).

Parametry ve funkci residuálu

Úlohu jsme vyřešili pro zadané složení kapaliny a zadaný tlak. Tyto parametry jsou napevno zapsany ve funkci residuálů res_boiling. Pokud se změní zadání, bude potřeba funkci residuálů přepsat. V této části si ukážeme, jak vytvořit funkci, která bude univerzálně použitelná pro různá složení kapaliny a různé tlaky.

Koncepce, kterou budeme často používat, bude založena na uložení všech parameterů modelu do jedné proměnné, tzv. vektoru parametrů, a předání tohoto vektoru z hlavního programu do funkce residuálu. Nejdříve si napišme seznam parametrů, které budeme potřebovat

tlak,
molární zlomek hexanu v kapalné fázi a
molární zlomek heptanu v kapalné fází.

Molární zlomek oktanu nemusíme zadávat, protože jej vždy můžeme dopočítat. Napišme teď upravenou verzi funkce res_boiling, která bude mít dva vstupní parametry: teplotu (stavovou proměnnou) a vektor parametrů

function res = res_boiling(t, pars)
%RES_BOILING  Residual function for boiling temperature calculation problem
%  Boiling temperature calculation by solving non-linear equation
%    t is temperature in Celsius
%    pars is vector of parameters
%      pars(1) - pressure in kPa
%      pars(2) - liquid mole fraction hexane
%      pars(3) - liquid mole fraction heptane
%    Raoult law for liquid/vapour equilibrium assumed
 
% unpack model parameters
patm = pars(1);    % atmospheric pressure in kPa
x1 = pars(2);
x2 = pars(3);
x3 = 1 - x1 - x2;  % liquid phase mole fractions
 
[A,B,C] = ptension_coef(); % get Antoine equation coefficients 
p1 = ptension(A(1), B(1), C(1), t); % hexane vapour pressure in kPa
p2 = ptension(A(2), B(2), C(2), t); % heptane vapour pressure in kPa
p3 = ptension(A(3), B(3), C(3), t); % octane vapour pressure in kPa
 
% calculate the residual
res = p1*x1 + p2*x2 + p3*x3 - patm;
end

Všimněte si, že jsme na začátek kódu přidali sekci “rozbalení parametrů”, kterou jsme nahradili dříve napevno zadané hodnoty.

Pokud budeme chtít použit odkaz na tuto funkci jako parametr při volání řešiče fzero, narazíme na problém. Funkce fzero vyžaduje odkaz na funkci pouze s jednou proměnnou, zatímco námi upravená funkce res_boiling má teď parametry dva. Toto je původní volání, které teď nefunguje

>> fzero(@res_boiling, 50)

Potřebujeme další funkci, která bude sloužit jako “adaptér” a umožní připojit res_boiling k fzero.
Bude jím anonymní funkce, která volá funkci res_boiling se známým vektorem parametrů a odkaz na tuto anonymní funkci pošleme jako první parametr funkci fzero

>> fzero(@(x) res_boiling(x,pars), 50)

Nakonec ukážeme, jak může vypadat hlavní program (skript), který řeší danou úlohu včetně výpočtu složení parní fáze.

% Main script for boiling temperature calculation problem
 
patm = 101.325; % pressure in kPa
x1   = 0.3333;  % hexane mole fraction in liquid
x2   = 0.3333;  % heptane mole fraction in liquid
x3 = 1-x1-x2;
 
% Pack vector of parameters
pars = [patm, x1, x2];
 
% initial guess for boiling temperature (Celsius)
t0 = 50; 
 
% Solve nonlinear equation for uknown boiling temperature
tb = fzero(@(x) res_boiling(x,pars), t0)
 
% Check that solution converged correctly
res = res_boiling(tb, pars)
 
% Calculate vapour pressure at boiling temperature
[A,B,C] = ptension_coef();
p1 = ptension(A(1),B(1),C(1),tb);
p2 = ptension(A(2),B(2),C(2),tb);
p3 = ptension(A(3),B(3),C(3),tb); % kPa
 
% Calculate vapour mole fractions
y1 = x1*p1/patm
y2 = x2*p2/patm
y3 = x3*p3/patm
 
% Check mole fraction sum equals one
sum_y = y1+y2+y3

Význam všech řádků skriptu by měl být jasný. Řešení soustavy více nelineárních rovnic si vysvětlíme příště.

Kontrolní otázky a úkoly k procvičení

Najděte teplotu varu a složení parní fáze pro kapalinu obsahující 12 mol% hexanu, 24 mol% heptanu a 64 mol% oktanu a tlak 5 bar.
Určete rosný bod a odpovídající složení kapalné fáze. Předpokládejte atmosférický tlak a ekvimolární složení parní fáze (molární zlomek hexanu, heptanu i oktanu je 1/3). Odvoďte potřebnou nelineární rovnici a napište k ní funkci residuálu res_dewpoint.m. Vyřešte nelineární rovnici, abyste nalezli teplotu rosného bodu a potom doočítejte složení kapaliny. Opakujte řešení pro jiné zadání. Přesvědčte se, že když zadáte složení páry stejné jako výsledek výpočtu bodu varu z příkladu 1, bude složení kapaliny odpovídat zadání příkladu 1.

Print This Post