Simulace – 6. díl

V elektrotechnické praxi se často objevuje napájení střídavým proudem. Jeho časový průběh je ne náhodou sinusový a my se dnes podíváme na jeden ze způsobů, jak jej vytvořit.

Funkce sinus je definována hned několikrát, podle stupně obtížnosti, na kterém se právě nacházíte. Na oblíbené Wikipedii lze najít definice z pythagorejského trojúhelníku, ale to vás jistě učili už na základní škole. Trošku složitější definici můžeme nalézt v Eulerově vzorci.

Nechci se tu zabývat nějakou složitou matematikou, jen si chci připravit půdu pro další výstavbu vědomostí. Vzorec totiž nápadně připomíná všeobecné řešení Cauchyho počáteční úlohy druhého řádu (se speciální pravou stranou) a právě toho využijeme při generování sinusového průběhu – budeme řešit speciálně upravenou diferenciální rovnici. Úprava je zde nutná, protože řešení systému (rovnice) druhého řádu obecně dává mnoho možných řešení, a tedy i různých průběhů odpovídajících vlastnostem systému: na stejnou otázku „Jak se máš?“ odpoví Pepa jinak než Láďa. Abychom se ale měli čeho chytit, budeme muset chvilku teoretizovat a rovnici si vyřešit obecně.

Zase ta matematika…

Náš systém tedy bude druhého řádu a jeho matematický model bude mít tvar

Druhý tvar obsahuje pouze zjednodušené označení derivací, nic víc. Tento systém vyřešíme následovně:

› Vypočítáme kořeny charakteristické rovnice.

› Zapíšeme všeobecné řešení:

Konstanty K1 a K2 vypovídají o amplitudě průběhu v přechodném ději a konstanta K3 odpovídá ustálenému stavu systému. Z následujícího zápisu se pokusíme určit vlastnosti systému. Očekáváme přechodný kmitavý děj, který neskončí, tj. ustálená složka v něm vystupovat nebude (K3 = 0).

Hledáme tedy takové kořeny charakteristické rovnice, pro které můžeme potvrdit předchozí rovnost. Hledáme vlastnosti systému, které vyvolají sinusovou odezvu na vnější podnět.

› Porovnáme obě strany rovnice:

Navíc víme, že kořeny charakteristické rovnice jsou komplexně sdružená čísla (i, −i). Pokud se podíváme na výpočet kořenů charakteristické rovnice, získáme komplexně sdružená čísla tehdy, když koeficient A bude nulový a zároveň koeficient B kladný.

Jen tehdy docílíme toho, že výraz pod odmocninou bude záporný a odmocninou ze záporného čísla je číslo komplexní. Koeficienty K1 a K2 získáme až dosazením pravé strany diferenciální rovnice. To ale pro nás není tak podstatné, získali jsme totiž tvar diferenciální rovnice, který na vnější podnět reaguje sinusovým průběhem jeho veličin (správně řečeno stavových veličin). Budeme tedy řešit rovnici ve tvaru:

Daná rovnice je druhého řádu, naším dalším úkolem tedy bude upravit rovnici tak, abychom ji dokázali vyřešit pomocí numerických metod, které jsou určeny k řešení rovnic prvního řádu. Na tuto úpravu se používá tzv. substituční (dosazovací) metoda. Zvolíme si substituci tak, aby v rovnici vystupovala derivace pouze prvního řádu. Samozřejmě, nic není zadarmo, takže naše substituce bude vystupovat jako další rovnice prvního řádu. Přesto jsme se ale dostali tam, kam potřebujeme – řešíme rovnice prvního řádu. Jsou sice dvě, ale to nám nevadí.

Máme tedy soustavu dvou rovnic prvního řádu. To, že se v první rovnici vyskytuje y(x) a ve druhé u(x), nás trápit nemusí. V každém kroku budeme tyto hodnoty znát z předchozího výpočtu. Nezbývá tedy už nic jiného, než že si rovnice přepíšeme tak, abychom je mohli přepsat do programu:

Skáče, jak já pískám

Jestliže jsme vytvořili samočinně kmitající soustavu, není potřebné ji budit vnějším podnětem, tedy C = 0. Vzpomeňme si na prosté kyvadlo: stačí, abychom jej jednou vychýlili, a kyvadlo začne samo kmitat. Pokud bychom neuvažovali tření o vzduch, kmitalo by kyvadlo pravděpodobně dodnes bez toho, abychom mu v tom pomáhali. Ve skutečném životě tento fakt zanedbat nemůžeme, ale v matematice si to dovolíme. Pokud si vzpomínáte, řešíme tzv. Cauchyho počáteční úlohu, tedy úlohu, kde známe stav na začátku a z něj určujeme stavy následující. Náš počáteční stav bude tedy jakýmsi vychýlením kyvadla, určením nenulové počáteční hodnoty stavových veličin.

Poslední zjednodušení se bude týkat koeficientu B. Pro první výpočty můžeme uvažovat, že B = 1. Tím bychom měli získat sinusový výstup s jednotkovou frekvencí, tedy jakýsi základ. Další úpravy pro změnu frekvence si vyzkoušíme až v hotovém programu.

Napíšeme tedy kroky, které musíme v každém cyklu výpočtu udělat:

Napíšeme si jednoduchý program v jazyku C. Záměnu indexů (poslední dva řádky předchozího výpisu) programovat nemusíme, protože v každém cyklu ukládáme výsledek do jedné proměnné a zároveň do souboru, proměnnou v dalším cyklu přepíšeme novou hodnotou.

// Generování funkce sinus do souboru sim.txt
// www.programujte.com 2006 - Matěj Pácha
// Visual C++ 2005 Express Edition

#include "stdafx.h"
#include 

using namespace std;

int _tmain(int argc, _TCHAR* argv[])
{
	ofstream fout;
	fout.open("sim.txt");

	double x,y,u,h,b,konec;
	x=0;			//nezávislá proměnná - odpovídá proměnné času
	y=1;			//nastavíme počáteční podmínky
	u=0;			//
	h=0.001;		//krok, se kterým počítáme
	konec=6.3;		//konec simulace (=2 pi)
	b=1;			//

	for(x=0; x<konec ; x+=h)
	{
		u += h*(-b*y);
		y += h*u;
		fout << x << "	" << u << "	" << y << "
";
	}
	fout.close();
	return 0;
}

Pokud budete chtít vykreslovat data ze souboru např. v tabulkovém editoru, budete pravděpodobně muset změnit v národním prostředí formát čísel (desetinná tečka). Na Internetu jsem ale našel zajímavý prográmek gnuplot.info, pomocí kterého jsem následující průběhy vykreslil bez změny formátu čísel:

Screenshot z gnuplotu:

Průběh veličiny u = f(x):

Průběh veličiny y = f(x):

Jak vidíte, průběh je sice sinusový, ale je časově (fázově) posunutý. Vše určuje počáteční podmínka. Zkuste si zvolit počáteční podmínku jinou, např. y(0) = −1 a uvidíte, že průběh veličiny u(x) bude sinusový s počátkem v nule. Uvedl jsem, že takto vykreslená funkce sinus má jednotkovou frekvenci. Tato frekvence však není časová, ale úhlová – periodou není jedna sekunda, ale obvod jednotkové kružnice − 2 pi = 6,28. Pokud změníme koeficient „b“, změníme časovou konstantu systému a soustava bude kmitat s jinou frekvencí. Nejrozumnější je si vše vyzkoušet.

Proč nepoužívat vestavěné funkce?

Kdo z vás kdysi programoval kalkulačku, určitě jste použili knihovnu spojitý děj. Při změně frekvence totiž přecházíte z jedné funkce na druhou v jednom okamžiku, kdy ale není zaručené, že obě funkce mají v tomto okamžiku stejnou hodnotu. Nastane tedy skok v průběhu, který se může ve složitějších simulacích projevit i rozkmitáním celého systému.

Na obrázku se mění frekvence z jednotkové (žlutě) na 2,5násobek (zeleně), výsledek je zvýrazněn červenou barvou. Pokud tedy potřebujete pouze vypočítat hodnotu funkce sinus, určitě vám stačí funkce z knihovny math (cmath). Výpočet se provádí pomocí Taylorova rozvoje, a je tedy i dostatečně přesný. Pokud ale chcete skokově měnit frekvenci bez toho, aby se narušila spojitost systému, budete muset sáhnout po simulačních metodách. Vlastnosti funkce sinus vypočítané simulací lze měnit v libovolném okamžiku. Jestliže chcete z nějakého důvodu generovat sinus fázově posunutý, stačí uvést obě počáteční podmínky nenulové. Pokud byste si vykreslili obě stavové veličiny u(x) a y(x) do jednoho grafu, zjistili byste, že průběhy jsou vzájemně posunuté o 90°. Tento posun platí vždy, a můžeme tedy podle něj určit počáteční podmínky pro sinus s libovolným fázovým posunutím. Amplitudu výstupního signálu můžeme měnit až ve výsledku pouhým vynásobením. Stejného efektu lze dosáhnout, pokud jako počáteční podmínky uvedeme u(0) = 0 a y(0) = max, kde max je žádaná amplituda signálu.

Ukázali jsme si způsob, jak lze generovat funkci sinus. Tímto způsobem se sinus generuje v případech, kdy hrozí ztráta stability systému. Já osobně jsem našel uplatnění v simulacích regulačních obvodů pro asynchronní motory, kdy se frekvence i amplituda mění na základě požadavku regulátoru. Funkcí sinus se ale také dá popsat pohyb bodu po kružnici, ačkoliv běžně se více v praxi využívají úhlově orientované souřadné systémy (cylindrický, sférický). Také jsme si ukázali příklad ideálního tzv. harmonického oscilátoru, jehož reálnými představiteli mohou být kyvadlo, soustava pružiny a závaží nebo elektrický LC oscilátor. Dalším důležitým bodem byla soustava druhého řádu, s níž budeme pracovat i příště.