Tartalomjegyzék

Analízis III (2017) - I. Vektoranalízis

file:   anal3_vekanal_1.m
author: Polcz Péter <ppolcz@gmail.com>

Created on 2017.08.04. Friday, 15:15:13

Az elkövetkezendőkben a Matlab Symbolic Math Toolbox (SMT) által kínált vektoranalízis (MathWorks) függvényeit használjuk: curl, divergence, gradient, jacobian, hessian, laplacian, stb...

További általan definiált vektoranalízishez kapcsolódó segédfüggvények (vekanal_) megtalálható a demonstrations GitHub repository-ban, a lib/matlab mappában:

% SMT examples
web(fullfile(docroot, 'symbolic/examples.html'))

% SMT vector analysis
web(fullfile(docroot, 'symbolic/vector-analysis.html'))

Divergencia, rotáció, gradiens, Jacobi mátrix számítása

Egy vektormező divergenciájának kiszámítása és ábrázolása

syms x y real
r = [x;y];

2D vektormező

F = [
    sin(x)
    cos(y)
    ];
pcz_pretty('F(x,y)', F)

Output:

F(x,y) =
/ sin(x) \
|        |
\ cos(y) /

Divergencia szimbolikus kiszámítása

divF = divergence(F,r);
pcz_pretty('div F(x,y)', divF)

Output:

div F(x,y) =
cos(x) - sin(y)

Szimbolikus skalárfüggvény numerizálása és ábrázolása

Mindenekelőtt létrehozok egy gridet, amelynek pontjait egyesével be kell helyettesíteni a szimbolikus függvénybe.

[xx,yy] = ndgrid(linspace(-4,4,30));

A behelyettesítés megoldható a subs paranccsal, ami viszont igen lassú. Ehelyett a szimbolikus függvényt anoním függvénnyé (function handle-vé) alakítom.

divF_fh = matlabFunction(divF)

Output:

divF_fh =
    @(x,y)cos(x)-sin(y)

Majd az anoním függvénybe behelyettesítem az értékeket. Ez lényegesen gyorsabb. A objektumok dimenzióira azonban nagyon vigyázni kell: xx és yy változók két dimenziós tömbbök. Érdemes őket 1 dimenziós tömbökké alakítani (jelen esetben 2D tömbökként is működne, de ha vektormezővel akarom ezt megtenni akkor a 2D tömbbök nagyon bekavarnak).

xx_1D = xx(:);
yy_1D = yy(:);

Most fog megtörténni a behelyettesítés:

divF_num_1D = divF_fh(xx_1D, yy_1D);

Visszaalakítás 2D-be

divF_num = reshape(divF_num_1D, size(xx));

Ennek egyszerűsítése végett létrehoztam egy vekanal_subsmesh függvényt, mely pontosan ezt a műveletet végzi el. Így kellene meghívni:

divF_num = vekanal_subsmesh(divF, r, {xx,yy});

figure('Position', [ 672 550 1149 430 ])
subplot(121), mesh(xx,yy,divF_num), hold on
title(sprintf('$F = %s, ~~ \\nabla F = %s$', latex(F'), latex(divF)), ...
    'interpreter', 'latex')
subplot(122), surf(xx,yy,divF_num, 'facealpha', 0.5), shading interp, grid off
view(0,90), axis equal, axis tight, hold on, colorbar
drawnow

Szimbolikus vektormező numerizálása és ábrázolása

Most ábrázoljuk a vektormezőt (majdnem ugyanaz a procedúra, egy kicsit nehezebb).

[F1,F2] = vekanal_subsmesh(F, r, {xx,yy});

subplot(121), quiver3(xx,yy,divF_num,F1,F2,F2*0,1.5)
subplot(122), quiver3(xx,yy,divF_num,F1,F2,F2*0,1.5)
drawnow

Deriváltak számítása

syms x y z real
r = [x;y;z];

% Vektormező
F = [
    x^2*y
    y*z
    x*y*z^2
    ];

% Skalármező
f = sin(x)*cos(y)*z^2;

pcz_display('F(x,y)', F)
pcz_display('div F(x,y)', divergence(F,r))
pcz_display('curl F(x,y)', curl(F,r))
pcz_display('jacobian F(x,y)', jacobian(F,r))

pcz_display('jacobian f(x,y) = `grad f(x,y)` [row-vector]', jacobian(f,r))
pcz_display('gradient f(x,y) [column-vector]', gradient(f,r))

Output:

F(x,y) [3x1] =

   x^2*y
     y*z
 x*y*z^2

div F(x,y) [sym] = z + 2*x*y + 2*x*y*z

curl F(x,y) [3x1] =

 x*z^2 - y
    -y*z^2
      -x^2

jacobian F(x,y) [3x3] =

[ 2*x*y,   x^2,       0]
[     0,     z,       y]
[ y*z^2, x*z^2, 2*x*y*z]

jacobian f(x,y) = `grad f(x,y)` [row-vector] [1x3] =

[ z^2*cos(x)*cos(y), -z^2*sin(x)*sin(y), 2*z*cos(y)*sin(x)]

gradient f(x,y) [column-vector] [3x1] =

  z^2*cos(x)*cos(y)
 -z^2*sin(x)*sin(y)
  2*z*cos(y)*sin(x)

Vektoranalízis deriválási szabályainak ellenőrzése (bonyolultabb)

Először deklaráljuk a szimbólikus objektumokat

syms x y z real
r = [x;y;z];

f = sym('f(x,y,z)');
g = sym('g(x,y,z)');

F = [
    sym('F1(x,y,z)')
    sym('F2(x,y,z)')
    sym('F3(x,y,z)')
    ];

G = [
    sym('G1(x,y,z)')
    sym('G2(x,y,z)')
    sym('G3(x,y,z)')
    ];

% Feltételezzük, hogy minden szimbólikus objektum valós (nem komplex)
assume([F,G,f,g,r],'real')

% Gradiens, divergencia, rotáció, vektoriális szorzat műveletek
grad = @(f) jacobian(f,r);
trans = @(F) reshape(F,[1 3]);

% grad: sorvektor
% gradient: oszlopvektor
% cross: olyan lesz a vegeredmeny mint amilyen az elso

% Végül a szabályok ellenőrzése, nullát kell kapjunk minden esetben
pcz_symeq(grad(f*g), f*grad(g) + g*grad(f), '∇(f g) = f ∇g + g ∇f')
pcz_symeq(divergence(f*F,r), grad(f)*F + f*divergence(F,r), '∇(f F) = ⟨∇f,F⟩ + ⟨f,∇F⟩')
pcz_symeq(curl(f*F), cross(gradient(f),F) + f*curl(F), '∇×(f F) = ∇f×F + f ∇×F')
pcz_symeq(divergence(cross(F,G),r), G'*curl(F) - F'*curl(G), '∇(F×G) = ⟨G,∇×F⟩ - ⟨F,∇×G⟩')
pcz_symeq(grad(F'*G), F'*grad(G) + G'*grad(F), '∇⟨F,G⟩ = ⟨F,∇G⟩ + ⟨G,∇F⟩')

Output:

[  OK  ] ∇(f g) = f ∇g + g ∇f
[  OK  ] ∇(f F) = ⟨∇f,F⟩ + ⟨f,∇F⟩
[  OK  ] ∇×(f F) = ∇f×F + f ∇×F
[  OK  ] ∇(F×G) = ⟨G,∇×F⟩ - ⟨F,∇×G⟩
[  OK  ] ∇⟨F,G⟩ = ⟨F,∇G⟩ + ⟨G,∇F⟩

Numerikus integrálás

Az ode45 függvény

Az ode45 egy numerikus solver, mely a $\frac{\mathrm{d} x}{\mathrm{d} t} = f(t,x)$ közönséges nemlineáris differenciál egyenletet oldja meg negyed- vagy ötödrendű Runge-Kutta módszerrel, változó időbeni $T_s$ lépésközzel. A differenciál egyenletben az ismeretlen egy $x: \mathbb{R} \mapsto \mathbb{R}^n$ időtől függő vektor értékű függvény.

Ezt a megoldót numerikus primitívfüggvény keresésére is használhatjuk. Legyen $F:\mathbb{R} \mapsto \mathbb{R}^n$ az ismeretlen vektor értékű $x$-ben változó függvény. Ismerjük $F$ függvény deriváltját: $F'(x) = f(x)$ és $F$ értékét $x_0$ pontban. Ekkor $F(x) = \int_{x_0}^{x}f(s)\mathrm{d}s$.

f = @(x,~) 2*x; % ez lehetne valami extémebb függvény is
x0 = 0; x_max = 5; F0 = -1;
[x,F] = ode45(f,[x0 x_max],F0);
figure, plot(x,F)

Az ode45 által visszaadott $m\times n$ méretű F tömb az $F(x)$ primitív függvény értékeit fogja tartalmazni az $x_i$, $i\in\{1,...,m\}$ pontokban, az $m\times 1$ méretű x tömb pedig az $x_i$ értékeket.

Az integral függvény

Az integral(fun,a,b) függvény kiszámítja egy |fun| anonim függvénnyel adott egyváltozós függvény határozott integrálját $D = [a,b]$ intervallumon. Számoljuk ki a következő integrált: $$\int\limits_0^{\frac{\pi}{2}} 4\cos(x)\sin(x)^2 - \cos(x)^2\sin(x) ~\mathrm{d}x$$

f = @(x) 4*cos(x).*sin(x).^2 - cos(x).^2.*sin(x);
a = 0;
b = pi/2;
I = integral(f, a,b)

Output:

I =
    1.0000

Vonalintegrál felületintegrál

Matlabban számítsuk ki numerikusan mekkora munkát gyakorol egy $\vec F(\vec r)$ erőtér egy anyagi pontra, miközben az anyagi pont a $\Gamma$ által leírt görbe mentén halad. Magyarán szólva, integráld az $\vec F(\vec r)$ vektormezőt $\Gamma$ görbe mentén. $\vec F(\vec r)$ és $\Gamma$ tetszőlegesen megválasztható, támpontként szolgálhat az 1. gyakorlat 5-7. feladata. Ugyanazokat a függvényeket is lehet használni.

Támpont

1. Definiáljunk egy $x$ és $y$-tól függő F, majd egy $t$-től függő gamma szimbolikus függvényt.
2. $\gamma(t)$ deriváltjának kiszámításához használjuk a diff parancsot.
3. A subs parancs segítségével helyettesítsük be [x;y] helyébe a gamma-t, hogy megkapjuk $\vec F(\vec \gamma(t))$-t.
4. A szimbolikus t változótól függő integrandust alakítsuk anoním függvénnyé a matlabFunction segítségével. Nevezzük el (pl.) Integrand_fh-nek
5. Végül hívjuk meg az integral parancsot numerikus integrálás végett:
   I = integral(Integrand_fh,t0,t1)

Matlabban számítsuk ki numerikusan mekkora a hozama egy $\vec F(\vec r)$ vektormező által leírt áramlásnak egy $\vec S(u,v)$ paraméteresen megadott felületen keresztül.
Vagyis: \begin{align} \int_S \vec F(\vec{r}) \mathrm{d}\vec{S} = \int_S F(\vec{r}) \vec{n}\mathrm{d} S = ? \end{align}

Támpont

Az eljárás majdnem ugyanaz mint vonalintegrál esetén, kettős integrál kiszámításához pedig használjátok az integral2 függvényt.

Vonalintegrál kiszámítása szimbolikus határozatlan integrálás segítségével

Integráljuk $\vec F(x,y)$-t $\Gamma$ mentén.

syms t x y a b real
r = [x;y];

$\vec F(x,y) = \vec F(\vec r)$ kétdimenziós vektormező

F = [
    x^2
    y^2
    ];

$\Gamma$ origó középpontú negyed ellepszis. Tengelyeinek hossza $2a$ és $2b$. Kezdőpont $A(a,0)$, végpont $B(0,b)$.

g = [
    a*cos(t)
    b*sin(t)
    ];

% Integrálási tartomány:
t1 = 0;
t2 = pi/2;

Integrálandó függvény: $\left<\vec F\big(\vec \gamma(t)\big),\dot{\vec \gamma}(t)\right>$

Integrand = subs(F,r,g)' * diff(g, t);
pcz_pretty('Integrálandó függvény', Integrand)

Output:

Integrálandó függvény =
 3              2    3       2
b  cos(t) sin(t)  - a  cos(t)  sin(t)

Szimbolikus határozatlan integrálás

Integral = int(Integrand, t)

Output:

Integral =
(b^3*sin(t))/3 + (a^3*cos(t)^3)/3 - (b^3*cos(t)^2*sin(t))/3

Határok behelyettesítése

I = subs(Integral,t,t2) - subs(Integral,t,t1)

Output:

I =
b^3/3 - a^3/3

Matlabos ügyeskedéssel így is lehet. FIGYELEM: Integral(t) nem egy anoním függvény. Integral(t1) esetén implicite meghívódik a subs parancs.

Integral(t) = int(Integrand, t)
I = Integral(t2) - Integral(t1)

Output:

Integral(t) =
(b^3*sin(t))/3 + (a^3*cos(t)^3)/3 - (b^3*cos(t)^2*sin(t))/3
I =
b^3/3 - a^3/3

Vonalintegrál kiszámítása szimbolikus határozott ill. numerikus módszer segítségével

syms t x y z real
r = [x;y;z];

$\vec F(x,y,z) = \vec F(\vec r)$ háromdimenziós vektormező

F = [
    3*x^2*y^2*z
    2*x^3*y*z
    x^3*y^2
    ];

$\Gamma$: Az origót az (1,2,3) ponttal összekötő szakasz.

g = t * [1;2;3];

% Integrálási tartomány:
t1 = 0;
t2 = 1;

Integrálandó függvény: $\left<\vec F\big(\vec \gamma(t)\big),\dot{\vec \gamma}(t)\right>$

Integrand = subs(F,r,g)' * diff(g, t)

Output:

Integrand =
72*t^5

Szimbolikus határozatlan integrálás

Integral = int(Integrand, t)

Output:

Integral =
12*t^6

Szimbolikus határozott integrálás

Integral = int(Integrand, t, t1, t2)

Output:

Integral =
12

Polinomközelítéses numerikus eljárással (először anoním függvénnyé kell alakítani)

Integrand_fh = matlabFunction(Integrand);
Integral = integral(Integrand_fh, t1, t2)

Output:

Integral =
    12

Felületintegrál számítása

syms u v real
syms x y z real
r = [x;y;z];

Vektormező

F = [
    x^2*sin(z)
    x*y
    z^2*x
    ];

Felület paraméteres megadása

S = [
    sin(u)*cos(v)
    sin(u)*sin(v)
    cos(u)
    ];

$\mathrm{d}\vec{S} = \left(\vec S'_u \times \vec S'_v \right) \mathrm{d}(u,v)$ számítása

dSu = jacobian(S,u)
dSv = jacobian(S,v)
dS = simplify(cross(dSu, dSv))

Output:

dSu =
 cos(u)*cos(v)
 cos(u)*sin(v)
       -sin(u)
dSv =
 -sin(u)*sin(v)
  cos(v)*sin(u)
              0
dS =
 cos(v)*sin(u)^2
 sin(u)^2*sin(v)
      sin(2*u)/2

$\vec F(\vec r)$-be helyettesítjük $\vec S(u,v)$-t.

FSuv = subs(F,r,S)

Output:

FSuv =
 sin(cos(u))*cos(v)^2*sin(u)^2
        cos(v)*sin(u)^2*sin(v)
        cos(u)^2*cos(v)*sin(u)

Itegrálandó: $\left<\vec F\big(\vec S(u,v)\big), \vec S'_u \times \vec S'_v\right>$

Integrand = FSuv' * dS

Output:

Integrand =
sin(cos(u))*cos(v)^3*sin(u)^4 + cos(v)*sin(u)^4*sin(v)^2 + (sin(2*u)*cos(u)^2*cos(v)*sin(u))/2

Anoním függvénnyé alakítás

Integrand_fh = matlabFunction(Integrand,'vars',{u v})

Output:

Integrand_fh =
    @(u,v)sin(cos(u)).*cos(v).^3.*sin(u).^4+cos(v).*sin(u).^4.*sin(v).^2+sin(u.*2.0).*cos(u).^2.*cos(v).*sin(u).*(1.0./2.0)

Integrálás $u = \left[\frac{\pi}{6},\frac{\pi}{3}\right]$, $v = \left[\frac{\pi}{4},\frac{\pi}{2}\right]$ tartományon.

integral2(Integrand_fh,pi/6,pi/3,pi/4,pi/2)

Output:

ans =
    0.0617

Published with MATLAB® R2015a, and sanitarized by BeautifulSoup (Python script using BSoup is written by Polcz).