[from Anthony] progress cont...

Just a quick comment on the asymptotic behaviour of D...what I thought for k3_real (red line) appeared as it looked was that we need to take care of the sign when doing a square-root...so one root for k=+sqrt(-D) and another for k=-sqrt(-D) which was why k3_real looked like a reflection of k1_real for r>>1. I agree...we need to do more analysis to explain k2 and also k3_real for r<<1.

****

So I tried integrating starting from the left with the following initial settings...

m=2 (yes...m was equal to 2)

dr = 0.003 (small enough to see if there could be kinks or interesting lumps)

r_initial = 0.01

eta(r=r_initial)= 0

v_eta(r=r_initial)=1

The coefficients B(r) and C(r) should be correct as I did the asymptotic predictions analytically and confirmed it numerically. Oh, from what I found, for r>>1, C(r) is actually increasing (rather than decreasing from our other parameters) because of the last term D/c^2, which is directly proportional to r (assuming c ~ v_kep). Nevertheless, the only part I think that has the greatest contribution to the solution of eta is I(r). Doing the asymptotic predictions for I(r) was a problem for me because I was not sure what assumptions I can do to simplify the term so I had to do it numerically.

Before I show the results, I'm going to define some of the symbols I used...

tau_1 ~ perturbed optical depth

c ~ sound speed

v_tau ~ dtau/dr (for unperturbed optical depth)

f ~ unperturbed radiative force

WW ~ shorthand for m*(Omega-Omega_p)

sigma ~ surface density

a ~ 2nd derivative of eta

v ~ 1st derivative of eta

From the list: sigma, c, f, WW, and v_tau have already been calculated prior to this integration.

So with the initial settings I adapt the leapfrog method for the integration...

tau_1[j+1] = tau_1[j]+(1/(c[j]*c[j]))*v_tau[j]*eta[j]*dr;

		

I[j+1]= I[j]+f[j]*(((log(r[j+1]*sigma[j+1]*f[j+1]/D[j+1])-log(r[j]*sigma[j]*f[j]/D[j]))/dr)+

				   ((2*m*W[j])/(r[j]*WW[j])))*tau_1[j];

		

a[j]= - B[j]*v[j] - C[j]*eta[j]- I[j];

		

v[j+1] = v[j] + a[j]*dr;

eta[j+1] = eta[j] + v[j+1]*dr;

So again, the only part I'm suspicious is tau_1 which includes I(r). But here's what I got from this....

I stopped the integration at R=5 because I thought it was meaningless to continue on. This seems to act very similar to the wave for k_3 (look back at my previous post). Now, I thought it was weird that 'nothing' was happening between r=0 to r=1 so I restricted r to this range and I get...

I noticed something was happening between r=0 and r=0.2 so I restricted r further...

Looking at the numerical results, it seems that C(r) was the dominating factor (which I think was the one making these responses) and after r=0.4, I(r) became the dominate one very quickly. Moreover, this also resembles the waves for k_1 and k_2. I also plotted eta between 0.4 to r=0.8-ish...

So if you piece together the last 3 pics you would get the first pic I posted.

That's all I have so far in terms of integrating. If you need to see more pics, let me know and I'll post them.

***

Pawel, I believe you reminded me that this second ODE we been dealing with is actually an integro-differential equation. So the only method that rings to my head to solve this kind is using Laplace transform. With the initial conditions we defined (eta(r_initial)=0; v_eta(r_initial)=1) it would be convenient to use this method where we transform and solve for L(eta) then inverse laplace transform it to get an analytic solution. The first couple of terms are easy to transform but for I(r) I was not sure whether we should treat (1/c^2)(dtau/dr) as a constant and pull them out of the integral. If I did that, it would be easier to solve. So now I'm at simplifying and solving for L(eta) and so far I have this fraction with a 3rd degree polynomial at the denominator. Now I just need to find some way to simplify it and inverse transform it. But before I go further, do you think this method is actually suitable for our case???

Anthony