Solving for a function, quite confounding

I am trying to solve for n_{H0} in the following physical equation:

Answer 1

[Roger Stafford]

Though you haven't said so, I assume that R, c, k, and T are also constants. I also assume that since a_nu as you defined it does not depend on the variable 'r', it can be removed from the inner integral inside the exponential, leaving just the integral of n_H0(r) with respect to r. Finally I assume that the r^2 in the denominator of the outer integral does not vary as nu is varied in the integration, so that it can be factored outside also. Are these all correct assumptions?

(By the way, it is poor practice to use the same symbol 'r' for the variable of integration on the inner integral and for the r^2 in the outer integral. Doing so can lead to confusion.)

Given these assumptions, this can be construed as a problem for matlab's 'fzero' function. Let us denote by 'x' the assumed value of the inner integral (without the a_nu) from R to infinity. From this assumed value of x we can do numerical integration with respect to nu to determine the outer integral (without the r^2). Then we can solve for n_H0 in the resulting quadratic equation and obtain it as a function of 'r' (which was removed from the outer integral.) Then numerical integration can be used on this function of r to obtain the actual inner integral of n_H0 with respect to r from R to inf which we could call y. The task for 'fzero' is to so adjust x that eventually y-x equals zero. There is of course the ambiguity that the quadratic equation has two roots, but very likely one of these can be quickly ruled out as physically inadmissible.

[Gloust]

Your first paragraph concerning assumptions are all correct.

As for your secondary comment, I should have written r prime(=r') so as to avoid possible confusion.

I am a tad confused about how to go about doing your description in the third paragraph. I understand we can take x:=integral from R to infinity of n_HO times dr. With this definition, I can do the outer integral to get some definite value. Then, if I read correctly, we solve this quadratic equation by the expression for quadratic roots(analytically). Solving by the quadratic root formula, I would then have an expression for n_H0 in terms of r and x. But I do not quite understand your explanation beginning with the sentence, 'Then numerical integration can be used on this function of r to obtain the actual inner integral of n_H0 with respect to r from R to inf...'

Thank you for your assistance, and hope you shall be able to clarify this for me.

[Roger Stafford]

"I understand we can take x:=integral from R to infinity of n_HO times dr. With this definition, I can do the outer integral to get some definite value." That isn't quite what I meant. I intended x to be some kind of guess on your part - or after the first cycle, on fzero's part - as to what this inner integral of n_H0 would turn out to be, not an actual integral evaluation. This would be the input to the function used by 'fzero'. First you would compute the outer integral based on this guess - you need the supposed value of the inner integral to solve the outer integral. Having done so, you would now be in a position to solve the resulting quadratic for n_H0 as a function of r (which was factored out of the outer integral.) Then finally you could use this function to compute the actual inner integral of n_H0 as a function of r, which you could call y. Naturally, it would in general be different from the earlier guess, x. The final value you would hand back to 'fzero' would be the difference y-x. If 'fzero' were being cooperative, there would come a time when the guess x becomes so accurate that it agrees with y as a result of your computation, at which time you have a solution, though conceivably there could be more than one such solution.

If I were doing your problem I might first let x range over some reasonable interval and compute y-x for each value of x. Actually you could compute two y's corresponding to each of the two quadratic roots. Then I would make a plot of both y-x's versus x to see at what value or values of x they appear to cross over zero. Then armed with this (or these) as initial estimates, I would put 'fzero' to work as outlined above. You should carefully read the documentation for 'fzero'.

I agree that this represents an unusually complicated procedure for computing the function given to 'fzero', but 'fzero' doesn't care how long it takes you to obtain the function value. It just blindly takes the output numbers it receives and attempts to adjust the input numbers so as to eventually get a zero output.

[Roger Stafford]

I have taken a closer look at your problem, Gloust, and if the interpretation I gave to your equation is correct, the equation cannot be satisfied unless the n_H constant is zero. The assumption that 1/r^2 factors out of the outer integral is what leads to that conclusion. The resulting quadratic equation would look like this:

   (n_H0-n_H)^2-I/alpha/r^2*n_H0 = 0

where I is the result of the outer integration and is presumably independent of r. It is obvious that as r approaches infinity, then the second term of the left side must approach zero and both roots of the quadratic equation in n_H0 must accordingly approach the constant n_H. If this latter is not zero, then the integrand of the inner integral which is also n_H0 will not be approaching zero as r approaches infinity, rendering that integral divergent.

In other words, either the assumptions made in that first paragraph that you stated were correct were actually incorrect in some respect, or else your problem has no solution unless n_H is equal to zero.

[Gloust]

I now noticed a slight mistake in the typesetting of the equation. The integral in the exponential is actually from 0 to r (instead of from R to infinity).

[Gloust]

n_H is constant, but as r approaches infinity, n_H0 approaches n_H (as there is the condition that n_H0+n_p=n_H, and n_p approaches zero as r approaches infinity). I did not include descriptors of all the physics in the exposition of the problem here, as it would have taken much space. Nevertheless, the equation is correct, I just am still very confused as to how to solve it explicitly via MATLAB for n_H0. I tried solving it but I am getting nonsensical distributions when I plot n_H0 vs. r. I am essentially getting n_H0 proportional to r^-2, which is unphysical because as we delve further from the star(go to bigger r), n_H0 should be increasing.

[Roger Stafford]

"The integral in the exponential is actually from 0 to r (instead of from R to infinity)." I'm afraid that still leaves you with a divergent inner integral. This time the difficulty occurs as r approaches the lower limit r = 0 in the inner integral. Because of the nature of your quadratic equation, the value of n_H0 is forced to approach infinity as r approaches 0, and that will make the integral divergent since it is essentially of a 1/r^2 nature which will give an infinite value to the integral.

I think you need to consider the possibility that your equation is not being interpreted correctly or possesses other fundamental errors.

[Gloust]

I would not be going quite to zero. I would start evaluating at say r=100, and this would be fine since the multiplicative factor of n_H0 in front(for which for small r is tremendously small) would suppress this tendency to approach big numbers. I still wonder if you may be so kind as to truly assist me with solving this numerically, as I tried it and as it is my first time trying such a thing my attempt my have been wrong.

[Gloust]

so the revised problem looks like this.