Solving a differential equation numerically to plot particle path

Question

I'm trying to plot the evolution of a particle in an accretion disk by solving the equation

$$2X\frac{\partial X}{\partial\tau}=V_R(X,\tau)$$

where I have found $V_R$ numerically to be

$$V_R(X_i,\tau_n)=-\frac{3}{2X_i^2\sigma_{i,n}}\bigg[\frac{X_{i+1}\sigma_{i+1,n}-X_{i-1}\sigma_{i-1,n}}{2\Delta X}\bigg]$$

and I have all values of the surface density $\sigma_{i,j}$.

So I assume I need to integrate $V_R$ with respect to time $\tau$ in order to get $X=f(\tau)$

I'm confused when it comes to trying to do this numerically. The first initial condition is $X(0)=\sqrt{0.9}$.

Any hints or help would be greatly appreciated because I've been stuck on this problem for way too long.

Answer 1

No, you don't want to integrate $V_R$ with respect to $\tau$. You want to solve the differential equation $2 X \frac{dX}{d\tau} = V_R(X,\tau)$ with initial condition $X(0) = \sqrt{0.9}$.

Unless you have very special requirements, you shouldn't try to re-invent the wheel. There is lots of good software for numerical solution of differential equations, and you can access it easily through Maple, Matlab, Mathematica etc.

Answer 2

To solve this numerically you would start by rearranging to get

$$\frac{\partial X}{\partial \tau} = \frac{V_R(X,\tau)}{2X} $$

and then starting with $X(0)$, you find $X(\Delta\tau)$ using:

$$ X(\Delta\tau) = X(0) + \frac{\partial X}{\partial \tau}\Delta\tau $$

and so on

$$ X((n+1)\Delta\tau) = X(n\Delta\tau) + \frac{\partial X}{\partial \tau}\Delta\tau $$

It is important that your $X$ variable be continuous. In other words, don't store it as $i$ or $X_i$, but rather as $X$ (otherwise you may find your particle doesn't move anywhere!). And then when you evaluate $\partial X/\partial\tau$ you can either round $X$ to the nearest $X_i$ or you can interpolate.