Changes between Version 2 and Version 3 of u/johannjc/injection


Ignore:
Timestamp:
12/15/17 15:49:48 (7 years ago)
Author:
Jonathan
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • u/johannjc/injection

    v2 v3  
    88First, let's consider a simpler example.  Let's assume we want to inject a certain amount of mass  $M_0$ around a point $\vec{X}$ (in 2D) with a kernel given by
    99
    10 $\delta \rho = M K(x)$
     10$\delta m = M K(x)$
    1111
    1212where our normalized kernel
     
    1818For a grid with 10 points in x and y over the interval x=[-1,1] y=[-1,1] - and for a sink particle at [.05,.05] the result is
    1919
    20 [[Image(Screen Shot 2017-12-15 at 3.06.47 PM.png​)]]
     20[[Image(Screen Shot 2017-12-15 at 3.10.01 PM.png​)]]
     21
     22Note the center of mass of the injected material is not .05,.05 and the injected mass is not 1.0.  This is because of discretization error.
     23
     24We would like to find a solution for $\delta m_i$ that is close to this, but subject to the constraints involving the center of mass and the total mass injected.  If we consider each zone inside the kernel, we can write the matrix equation
     25
     26{{{
     27#!latex
     28 \[
     29 \left[\begin{array}{rrrrr}
     30  x_1 & x_2 & x_3 & ... & x_n  \\
     31   y_1 & y_2 & y_3 & ... & y_n  \\
     32   1 & 1 & 1 & ... & 1  \\
     33   \end{array} \right]
     34 \left[\begin{array}{r}
     35  \delta m_1  \\
     36  \delta m_2  \\
     37  \delta m_3  \\
     38  ...  \\
     39  \delta m_n  \\
     40   \end{array} \right]=
     41 \left[\begin{array}{r}
     42  MX_0  \\
     43  MY_0  \\
     44  M  \\
     45   \end{array} \right]
     46\]
     47}}}
    2148
    2249
     50This can be expressed as $Ax=b$ and we can solve this underconstrained problem using the pseudo-inverse of A
    2351
    2452[[Image(Screen Shot 2017-12-15 at 3.07.00 PM.png​)]]
    2553
     54While this solves our constraints - it looks nothing like our desired profile.  What we would like to do is to solve $Ax=b$, in a way that does not minimize $||x||$, but that minimizes $||x-d||$ where $d$ is our vector of initial $\delta_m$ given by our kernel.
     55
     56We can do this using Matlab's lsqlin function - or we can use the prescription given in https://see.stanford.edu/materials/lsoeldsee263/08-min-norm.pdf (page 8-13 through 8-15)
     57
     58This gives
     59$x=d-A^T(AA^T)^{-1}(A d - b)$
     60
    2661[[Image(Screen Shot 2017-12-15 at 3.07.04 PM.png)]]
     62
     63
     64Now for the more complicated case, we have $\delta m_i, v_i$, but the matrix would contain two blocks - so we can solve the problem in two steps.  First we can solve for $\delta m_i$, and then use that solution to construct the matrix equation for $v_i$
     65
     66{{{
     67#!latex
     68 \[
     69 \left[\begin{array}{rrrrrrrrr}
     70  m_1 & ... & m_n & 0 & ... & 0 & 0 & ... & 0 \\
     71  0 & ... & 0 & m_1 & ... & m_n & 0 & ... & 0 \\
     72  0 & ... & 0 & 0 & ... & 0 & m_1 & ... & m_n  \\
     73   0 & ... & 0  & m_1z_1 & ... & m_nz_n &  -m_1y_1 & ... & -m_ny_n\\
     74  -m_1z_1 & ... & -m_nz_n &  0 & ... & 0 & m_1x_1 & ... & m_nx_n \\
     75m_1y_1 & ... & m_ny_n &  -m_1x_1 & ... & -m_nx_n &0 & ... & 0  \\ 
     76   \end{array} \right]
     77 \left[\begin{array}{r}
     78  v_{x1}  \\
     79  ...  \\
     80  v_{xn} \\
     81  v_{y1}  \\
     82  ...  \\
     83  v_{yn} \\
     84  v_{z1}  \\
     85  ...  \\
     86  v_{zn} \\
     87    \end{array} \right]=
     88 \left[\begin{array}{r}
     89  0  \\
     90  0  \\
     91  0  \\
     92  J_x  \\
     93  J_y  \\
     94  J_z  \\
     95   \end{array} \right]
     96\]
     97}}}