1
$\begingroup$

I am trying to solve a set of couple ODEs:

$V_l(r) - r W_l(r) - f1(r) W_l' = 0\tag 1$

$r^2 h''_l(r) + f2 r h_l'(r) + f3 h_l(r) - f4 U_l(r) = 0 \tag 2$

$\kappa (U_l + h_l) + V_{l+1} + W_{l+1} = 0\tag 3$

$\kappa (U_{l-2} + h_{l-2}) + f5V_l' + f6W_l' + f7 U'_{l-1} +f8U'_{l-2} = 0 \tag 4$

for $V_l,W_l,h_l,U_l$ and $\kappa$. These functions and the $fn$ are only functions of $r$. Note: These are not the actual ODEs, just a representation of them. There are actually not 4 equations, but $4*l$, where $l$ is the one that appears in the spherical harmonics $Y^l_m$, and $l$ must start at $l=2$. Also I know what the $fn$ functions are, I do not need to solve for them.

The method I have used to solve this is a spectral approximation. I take know and unknown functions in the ODEs and expand them in Chebyshev polynomials, for example

$f1(r) = \Sigma_{i=0}^{\infty}f1_i T_i[y] - \frac{1}{2} f1_0 $

and

$U_2 = \Sigma_{i=0}^{\infty}U_{2_i} T_i[y] - \frac{1}{2} U_{2_0} $

Applying this to equation (2) I get something like:

$[\kappa * (U_{2_N} + h_{2_N}) + V_{3_N} + W_{3_N} ]*T_N + \ldots - 1/2[\kappa * (U_{2_0} + h_{2_0}) + V_{3_0} + W_{3_0} ]*T_0 $

where $N$ is the maximum order of the Chebyshev polynomial. So if I extract the coefficients of the Chebyshev polynomials I get $4*l*N$, non-linear equations. I say nonlinear because there will be terms like $\kappa * U_{2_N}$, etc. From now on I will refer to the vector of $U_{2_N}, h_{2_N}, \ldots$ as $X$.

In order to solve this non-linear set of equations I am first using the condition that the matrix of equations, $a(\kappa,X)$ needs to satisfy $det(a) = 0$. To makes use of this I am using a singular value decomposition (SVD) to find a $\kappa$ that makes $a$ as singular as possible. To get the values for $X$ I make a new matrix $a'$ in which I remove a row from $a$ and replace it with a "normalization" condition

$\Sigma U_{2_i} T_i = 1$

and then do LinearSolve[a',B]. All the elements of B are zero, expect the one that corresponds to the row that I delete, that element of B is 1.

The main issue I am having is that my value for $\kappa$ does not converge very well. That is, when I increase $l$ and increase $N$ I would expect to get one value of $\kappa$, however at best I can get about 4 digits of convergence, but no matter what I try to do the fifth digit is just random

I've tried messing with the physical model, and trying to increase the accuracy of the calculations, but nothing seems to help. The matrix is poorly conditioned: 1^8 to 1^10. And so I'm sure this is causing some issues.

My attempts at a solution 1st Attempt:

I've tried to implement a Newton-Raphson method to try and fix the convergence. What I do is start with the matrix $a$, and treat $\kappa$ as a variable like one of the variables in $X$. Then I create the Jacobian

$J = \frac{\partial A_i }{\partial x_j}$ where $A_i$ represents an element of a vector that is made from the rows of my matrix $a$ with the additional equation that comes from the "normalization" condition. I include this so that I have $4*l*N +1 $ equations since now I'm trying so solve for $\kappa$ and $X$ all at once. And $x_j$ has all the elements of $X$ and $\kappa$.

Then I use the values of $\kappa$ and $X$ I got from doing the SVD and row replacement stuff above to use as the initial guess for the Newton's method, and then do

$z = LinearSolve[J, -A]$, and I take $z$ and add it to guess $X$, and repeat. However there is no convergence, in fact the the vector $X$ goes to zero like I'm stuck at the root, basically the $z$ I get is just $-X$. I tried giving LinearSolve the Krylov option, and there the $z$ that I get just flies off to infinity.

I'm not sure where the problem is. I'm thinking it has to do with the fact that my Jacobian might be close to singular ($det (J) \approx.00001$). I read some other stackexchange forums that said to use a QR decomposition instead of LinearSolve, and I tried that and got no convergence there either.

2nd Attempt: I tried using a QR algorithm for finding eigenvalues on matrix $a(\kappa,X)$. The idea is that I will look for the value of $\kappa$ which give me an eigenvalue as close to zero as possible. Again, this method does not do any better than the SVD stuff I was already using.

Long, long, story short, I really need some help, figuring out just how to make this result converge better. I would want something like 10 digits, but seriously 6 or 7 would be a miracle. What's more I don't care how slow the code gets, if it works it works.

1st Edit The functions are defined on the domain ($0,R$), where physically zero is the center of a sphere and $R$ is the radius.

The boundary conditions are: $W_l(r = R) = 0$ $h_{in}(r=R) = h_{out}(r=R)$

and $h_{l_{out}(r)} = constant*(R/r)^l$

and $h_l \rightarrow (r/R)^{l+1}$ as $r \rightarrow 0$

$U_l \rightarrow (r/R)^{l+1}$ as $r \rightarrow 0$

$V_l \rightarrow (r/R)^l$ as $r \rightarrow 0$

$W_l \rightarrow (r/R)^l$ as $r \rightarrow 0$

$\endgroup$
  • $\begingroup$ Please specify the boundary/initial conditions, so that we can tell whether you have an initial or boundary value problem (and what the domain is). Depending on which it is, there are plenty of standard methods available. $\endgroup$ – David Ketcheson Oct 26 '13 at 6:53
  • $\begingroup$ Thanks for looking at this David, I've adde the boundary conditions. $\endgroup$ – tau1777 Oct 26 '13 at 16:30
  • $\begingroup$ I see now that you have a boundary value problem. Have you tried just handing your problem to a standard BVP solver in MATLAB, Mathematica, or Maple? $\endgroup$ – David Ketcheson Oct 27 '13 at 12:05
  • $\begingroup$ I don't think I've ever really tried that, because of the nature of the boundary conditions. Not really sure how to give Mathematica something like $h''_l + h'_l +\ldots = 0$ given $h_l(r=0) \rightarrow (r/R)^3$. Also there is the issue of the parameter $\kappa$. I would have to guess a value and solve the ODEs and somehow keep guessing until things converged. Not sure how to do this either? Thanks again for looking at this. $\endgroup$ – tau1777 Oct 28 '13 at 20:30
  • $\begingroup$ I don't understand what you mean by $h_l(r=0)\to(r/R)^3$. Shouldn't that just be $h_l(r=0)\to0$? And I don't understand where the other condition comes from in your description. Anyway, I suggest re-writing your question with the whole problem statement at the top and an explanation of why you haven't used a standard method. $\endgroup$ – David Ketcheson Oct 29 '13 at 3:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.