Uvažujme rovnici $$ u_t + f(u)_x = 0 $$
s nespojitým (slabým) řešením ve tvaru $$ u(x,t) = \left\{ \begin{array}{ll} u_l & \text{pro\ } x < \xi(t), \\ u_r & \text{pro\ } x \ge \xi(t), \end{array} \right. $$ kde funkce $\xi(t)$ udává polohu nespojitosti v čase $t$. Potom pro rychlost pohybu nespojitosti $s = \dot{\xi}(t)$ platí $$ s = \frac{f(u_r)-f(u_l)}{u_r-u_l}. $$
In [1]:
using PyPlot
In [2]:
Nx = 101; # Pocet bodu site na intervalu [0,1]
x = range(0, stop=1,length=Nx); # Souradnice bodu site
function u0(x)
"""Definice pocatecni podminky"""
if x<=0.25
return x
else
return 0
end
end
plot(x,[u0(xi) for xi in x]);
ylim(-0.1,0.3); grid(true); xlabel("x"); ylabel("u");
In [3]:
dx = 1.0 / (Nx-1); # Velikost prostoroveho
dt = 0.5 * dx / 1.0 # Velikost casoveho kroku
Tend = 10;
pocet_iteraci = round(Int, Tend / dt) # Vypocet provedeme do casu priblizne Tend
t = zeros( pocet_iteraci )
u = zeros( (Nx, pocet_iteraci) )
u[:,1] = [u0(xi) for xi in x]
function flux_upwind(ul, ur)
if (ul+ur>0)
return ul^2/2
elseif (ul+ur<0)
return ur^2/2
else
return 0;
end
end
for n in 1:pocet_iteraci-1
for i in 2:Nx-1
u[i,n+1] = u[i,n] - dt / dx * ( flux_upwind(u[i,n],u[i+1,n]) - flux_upwind(u[i-1,n],u[i,n]) )
end
t[n+1] = t[n] + dt
end
plot(x, [u0(xi) for xi in x], label="u0")
T=1/4*Tend; plot(x,u[:,div((pocet_iteraci-1), 4)], "-", label="t=$T");
T=1/2*Tend; plot(x,u[:,div((pocet_iteraci-1) * 2,4) ], "-", label="t=$T");
T=3/4*Tend; plot(x,u[:,div((pocet_iteraci-1) * 3,4) ], "-", label="t=$T");
T=Tend; plot(x,u[:,pocet_iteraci-1], "-", label="t=$T");
ylim(-0.1,0.3); grid(true); xlabel("x"); ylabel("u"); legend(loc="upper right");
In [4]:
contourf(collect(x'), t, collect(u'));
xlabel("x"); ylabel("t");
In [5]:
#plot_surface(x', t, u', rstride=20, cstride=1, linewidth=0, cmap=ColorMap("rainbow"));
#xlabel("x"); ylabel("t"); zlabel("u");
Nechť $z \in \mathbb C^2(\mathbb R \to \mathbb R)$ je konvexní ($z''(u) \ge 0$) a $g'(u) = z'(u) f'(u)$. Potom pro hladké řešení rovnice $$ u_t + f(u)_x = 0 $$ platí $$ z(u)_t + g(u)_x = 0. $$
Pokud je (slabé) řešení nespojité, tak nemusí být jednoznačné. Jednoznačné řešení lze pak získat jako $$ u(x,t) = \lim_{\epsilon\to 0+} u^\epsilon(x,t), $$ kde $u^\epsilon$ je řešením rovnice $$ u^\epsilon_t + f(u^\epsilon)_x = \epsilon u^\epsilon_{xx}. $$
Pro takto získané řešení platí (ve slabém smyslu) $$ z(u)_t + g(u)_x \le 0. $$
Uvažujme Burgersovu rovnici $u_t + (\frac{1}{2}u^2)_x=0$ se slabým řešením $u(x,t) = 1$ pro $x<0$ a $-1$ pro $x \ge 0$. Zvolme $z(u) = u^2/2$. Poznámka: správně bychom měli entropickou podmínku ověřovat pro všechny konvexni funkce $z$, to lze nahradit vyšetřením pro speciální entropie $z^p(u) = max(u-p,0)$, kde $p\in\mathbb R$. My se však pro demonstraci omezíme na výše uvedenou entropii. Pro tuto entropii je $$ g'(u) = z'(u) f'(u) = u u = u^2 $$ a tedy $$ g(u) = \frac{1}{3} u^3. $$
Integrujeme-li entropickou nerovnost na intervalech $x \in [-1,1]$ a $t \in [0,1]$, dostáváme $$ \int_0^1 \int_{-1}^1 (z(u)_t + g(u)_x)\,dx\,dt = \int_{-1}^1 [z(u(x,t))]_{t=0}^1 \,dx + \int_{0}^1 [g(u(x,t))]_{x=-1}^1 \,dt = (*) $$ Vzhledem k tomu, že řešení je po částech konstantní, zjednoduší se nám výše uvedený integrál na $$ (*) = 2 (\frac{1}{2}-\frac{1}{2}) + 1 \frac{1}{3}( (-1)^3 - 1^3 ) = - \frac{2}{3} \le 0. $$
Uvažujme Burgersovu rovnici $u_t + (\frac{1}{2}u^2)_x=0$ se slabým řešením $u(x,t) = -1$ pro $x<0$ a $1$ pro $x \ge 0$. Postupujeme stejně jako v předchozím příkladě se $z(u) = u^2/2$.
Integrujeme-li entropickou nerovnost na intervalech $x \in [-1,1]$ a $t \in [0,1]$, dostáváme $$ \int_0^1 \int_{-1}^1 (z(u)_t + g(u)_x)\,dx\,dt = \int_{-1}^1 [z(u(x,t))]_{t=0}^1 \,dx + \int_{0}^1 [g(u(x,t))]_{x=-1}^1 \,dt = (*) $$ Tentokrát je však $$ (*) = 2 (\frac{1}{2}-\frac{1}{2}) + 1 \frac{1}{3}( 1^3 - (-1)^3 ) = \frac{2}{3} > 0. $$
In [6]:
Nx = 51; # Pocet bodu site na intervalu [0,1]
x = range(0, stop=1, length=Nx); # Souradnice bodu site
function u0(x)
"""Definice pocatecni podminky"""
if x<0.5
return -1.0
else
return 1.0
end
end
function uex(x,t)
"""Presne reseni"""
if (abs(x-0.5)<t)
return (x-0.5)/t
elseif(x<0.5)
return -1.0
else
return 1.0
end
end
plot(x,[u0(xi) for xi in x], label="u0");
plot(x,[uex(xi, 0.25) for xi in x], label="uex");
legend(loc="upper left");
ylim(-1.5,1.5); grid(true); xlabel("x"); ylabel("u");
In [7]:
dx = 1.0 / (Nx-1); # Velikost prostoroveho
dt = 0.5 * dx / 1.0 # Velikost casoveho kroku
u = [u0(xi) for xi in x]
uNew = copy(u)
t = 0
Tend = 0.25;
pocet_iteraci = round(Int, Tend / dt) # Vypocet provedeme do casu priblizne Tend
f(u) = u^2/2
function flux_upwind_1(ul, ur)
if (ul+ur>0)
return f(ul)
elseif (ul+ur<=0)
return f(ur)
end
end
for n in 1:pocet_iteraci
for i in 2:Nx-1
uNew[i] = u[i] - dt / dx * ( flux_upwind_1(u[i],u[i+1]) - flux_upwind_1(u[i-1],u[i]) )
end
t = t + dt
u = copy(uNew)
end
u_upwind = copy(u)
plot(x, [u0(xi) for xi in x], label="u0")
plot(x,[uex(xi, 0.25) for xi in x], label="uex");
plot(x,u, "o", label="u[i]");
ylim(-1.5,1.5); grid(true); xlabel("x"); ylabel("u"); legend(loc="upper left");
In [8]:
u = [u0(xi) for xi in x]
uNew = copy(u)
function flux_upwind_2(ul, ur)
if (ul+ur>0)
return f(ul)
elseif (ul+ur<0)
return f(ur)
else
return 0
end
end
for n in 1:pocet_iteraci
for i in 2:Nx-1
uNew[i] = u[i] - dt / dx * ( flux_upwind_2(u[i],u[i+1]) - flux_upwind_2(u[i-1],u[i]) )
end
t = t + dt
u = copy(uNew)
end
plot(x, [u0(xi) for xi in x], label="u0")
plot(x,[uex(xi, 0.25) for xi in x], label="uex");
plot(x,u, "o", label="u[i]");
ylim(-1.5,1.5); grid(true); xlabel("x"); ylabel("u"); legend(loc="upper left");
In [12]:
u = [u0(xi) for xi in x]
uNew = copy(u)
function myabs(s)
eps = 0.2
a = abs(s)
if (a>eps)
return a
else
return (a^2+eps^2)/(2*eps)
end
end
function flux_upwind_e(ul, ur)
if ul > ur
return flux_upwind_1(ul,ur)
else
s = (ul+ur) / 2.0
flux = f((ul+ur)/2.) - myabs(s)/2. * (ur-ul)
#flux = (f(ul) + f(ur))/2 - 0.5 * dx/dt/2 * (ur-ul)
return flux
end
end
for n in 1:pocet_iteraci
for i in 2:Nx-1
uNew[i] = u[i] - dt / dx * ( flux_upwind_e(u[i],u[i+1]) - flux_upwind_e(u[i-1],u[i]) )
end
t = t + dt
u = copy(uNew)
end
plot(x, [u0(xi) for xi in x], label="u0")
plot(x,[uex(xi, 0.25) for xi in x], label="uex");
plot(x,u, "o", label="u[i]");
ylim(-1.5,1.5); grid(true); xlabel("x"); ylabel("u"); legend(loc="upper left");
In [ ]:
In [ ]: