# Given a vector field construct a corresponding stream function

this might be a little complicated, so let's break it down to bits.

What I want to achieve

Firstly, every piece of data I'm referring to is numerical.

Let's presume we have a div-free (up to some prescribed precision) vector field $V = (u, v)$ (for simplicity, in two dimensions). This condition is crucial. Then there is a scalar (in 3D this would be also a vector) field $\psi$ defined as

\begin{aligned} u &= \frac{\partial \psi}{\partial y} \\ v &= - \frac{\partial \psi}{\partial x} \end{aligned}

Provided the vector field $(u,v)$ is div-free the formula for $\psi$ is: $$\psi (x, y) = \int \limits_{(x_0, y_0)}^{(x, y)} \left( u \, \mathrm{d} y - v \, \mathrm{d} x \right)$$

This expression only makes sense (not dependent on path taken in the integration) if $\partial_x u + \partial_y v = 0$.

I want Mathematica to calculate (numerically) stream function $\psi$ from given div-free vector field in two dimensions. I want to be able to choose the value $\psi_0 \equiv \psi (x_0, y_0)$.

Why I need it

Currently, I use ListStreamPlot for visualisation of vector fields. However, my vector field is changing in time and MMA's automatic choice of streamline positions (seed points) is changing abruptly over a time step. Maybe a better choice would be to seed the points for streamlines by hand, but I don't know where are the "interesting" parts of the vector field in advance! Moreover, streamlines produced by MMA are quite often ended in the middle of nowhere (I assume MMA is struggling with fields with a small magnitude). Streamlines (integral curves of a div-free vector field) are a solution to the equation

$$\psi (x, y) = const.$$

I presume that if I fixed the value of $\psi$ somewhere on the boundary, the contours (visualised by the ListContourPlot with some prescribed, constant Contours option) of this field would reflect only real changes in the corresponding vector field and not some arbitrary choices that MMA makes when trying to find interesting parts of the vector field. The Contours option would be chosen in a way the interesting parts would be visible.

My try on this

Currently, I am using the following code:

f[x_, y_] := x^3 - y^4 x;
v = {};
Monitor[For[i = 1, i <= 10000, i++,
With[{x = Random[], y = Random[], vx = D[f[r, s], s],
vy = -D[f[r, s], r]},
AppendTo[v, {{x, y}, {vx, vy} /. {r -> x, s -> y}}];
]
], i]
ListStreamPlot@v
ContourPlot[f[x, y], {x, 0, 1}, {y, 0, 1},
Contours -> Table[Sign[i] Abs[i]^1.6, {i, -1, 1, 0.01}],
WorkingPrecision -> 100, PlotPoints -> 100]

pts = Table[v[[i, 1]], {i, 1, Length@v}];
vx = Table[{v[[i, 1, 1]], v[[i, 1, 2]], v[[i, 2, 1]]}, {i, 1,
Length@v}]; vy =
Table[{v[[i, 1, 1]], v[[i, 1, 2]], v[[i, 2, 2]]}, {i, 1, Length@v}];
vx = Interpolation[vx];
vy = Interpolation[vy];
\[Psi] = {};
ClearAll[x0, y0, x, y];
Monitor[For[i = 1, i <= Length@pts, i++,
With[{x0 = 0, y0 = 0, x = pts[[i, 1]], y = pts[[i, 2]]},
val = NIntegrate[
vx[x0 + t (x - x0), y0 + t (y - y0)] (y - y0) -
vy[x0 + t (x - x0), y0 + t (y - y0)] (x - x0), {t, 0, 1}];
AppendTo[\[Psi], {x, y, val}];
]
], ToString@i <> "/" <> ToString@Length@pts]
ListContourPlot[\[Psi],
Contours -> Table[Sign[i] Abs[i]^1.6, {i, -1, 1, 0.01}]]


where $v$ is a list of pairs in the form: $v = \{ \{ \{ x_1, y_1 \}, \{ v_{x1}, v_{y1} \} \}, \{ \{ x_2, y_2 \}, \{ v_{x2}, v_{y2} \} \}, \cdots \}$.

Just for the record, the path integral mentioned above is implemented for linear path between points $(x_0, y_0)$ and $(x, y)$ interpolated by the variable $t$ changing from $0$ to $1$.

The problem is pretty obvious: firstly, the numerical integration is very slow. It takes several minutes to run the second part of the code (integration) on such a small sample (only 10k of points). The second one is precision:

I provided an example with exactly zero divergence, however, my fields (obtained numerically) always have although small, non-zero divergence. We can see the effect of this on the farthest points (from x0, y0). Apart from that, the results are pretty close.

Of course, another approach would be solving a poisson equation of the type

$$\Delta \psi = - \omega$$

where

$$\omega = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}$$

however, there are two problems:

The first problem is that my field $(u, v)$ is defined on an unstructured grid - there is only linear interpolation available in this case and I am worried about precision.

The second problem are boundary conditions for $\psi$ - what BC's should I apply to get the correct result? For this reasons I resorted to the path integral solution.

What would you suggest? What is the best approach to tackle this problem?