Solving a system of coupled iterative equations containing lots of heavy integrals and derivatives
Hello all,
I am trying to solve a system of coupled iterative equations, each of which containing lots of integrations and derivatives.
First I used Sage to solve it analytically, but the solution was too dependent on the initial guesses for my unknown functions in the iterative loops, constant values for initial guesses yielded in Sage answering back almost immediately whereas symbolic functions when were used as initial guesses caused the system to go deep into calculations, sometimes seemingly into never ending ones.
However, what I tried with Sage was actually a simplified version of my original equations so I thought it might be the case that I have no other choice rather than to treat the integrations and derivatives numerically, however, I had some not ignorable problems:
integrations were only allowed to have numerical limits and variables were not allowed as e.g. their upper limits (I thought maybe a numerical method algorithm is faster than the analytic one even-though I leave a variable or parameter in its calculations, but it just didn't work so).
integrands couldn't also admit extra variables and parameters w.r.t. which not being integrated.
the derivative function was itself a big obstacle, as I wasn't able to compute partial derivatives or to use a derivative in the integrand of an integral.
To get rid of all the problems with numerical derivative I substitute it by the symbolic diff() function and the speed improvement was still hopeful but the problems with numerical integration persist.
Now I have three questions:
a- Is it right to conclude there is no other way for me rather than to discretize the equations and do a complete numerical treatment instead of a mixed one?
b- If so then is there any way to do this automatically? My equations are not DE ones to use ODEint or else, they are iterative equations, I have integrations and derivatives only to update my unknowns at each step to their newer values.
c- If my calculations are so huge in size is there any suggestion on switching from python to fortran or things like that as well?
Best Regards
I have to solve big systems of coupled integro-differential equations for my work. I tried Mathematica and sage both are too slow for me. I ended up writing the code in c using gsl. If you are looking for a balance of speed and least work a c code is using gsl is the best option. I also tried using gsl routines in sage and running it in cython, but from my personal experience packages like sage are not meant for heavy duty specialized codes like the ones you are talking about. c or fortran is the best option. Are your equations linear or non-linear and how big is your system? I have solved upto 10 k coupled non-linear differential equations numerically using gsl in c. However, beyond that even gsl fails and I am in the process of writing my own algorithm for bigger systems.
thanks for your kind comment, I was indeed dealing with a huge system of coupled nonlinear integro-PDEs, however I was told it is not solvable using the present computers, so first I reduced them semi-analytically to some systems of iterative equations only containing plenty of integrations and derivatives, but even doing so Sage proved to be too slow. Here is my analytic try with Sage for my simplest case (12 equations and 12 unknowns): http://ask.sagemath.org/question/1661/how-to-speed-up-a-code-containing-several-symbolic , however things should get really more subtle when i reach the biggest case which is almost 4.5E(6) equations and unknowns. If this is not something that even gsl cannot handle then how should I develop my own algorithm?
Just by looking at your equations it seems like a very difficult problem. Do you know anything about what the solutions are supposed to look like? Are there solution that are rapidly oscillating as a result of turbulence that you are not interested in? It seems like that to me. If that is the case you will need to discretise it and you will get a lot of non-linear differential equations with the integrals replaced by sum. You can express it as df/dt= Af where f is a vector and A is a matrix. Then the best option for solving this would be the Magnus method (google it!). Basically you will find f = T(exp(At)) where T implies time ordering. In order to calculate exp(At) you will need eigenvalues and eigenvectors of A. If they are not calculable (numerically) then I don't know.
The solutions I expect are regular functions far from being rapidly oscillating in the sense of turbulence. Actually, I treat turbulence as a stochastic process, then expand its variables in a random space, the expansion coefficients of which then become my new unknowns which are by themselves deterministic, hopefully regular and smooth functions. I am not sure if they would be periodic at all for all available boundary conditions. I had hope to handle this (semi)analytically as it is what is really intended from semi-analytic methods like Adomian's decomposition method, but if there is no other way then I have no other choice. What about mixed analytical-numerical methods, e.g. using the numerical algorithms with some variables and parameters unassigned therein?
You will have a much better change if you know the eigenvalues and eigenvectors of the matrix 'A' I was talking about analytically or at least a part of it so that you can treat the rest as perturbation. That may be possible if the A is more or less sparse or block diagonal. The problem in your case is that the system is in 3D. So even if you make a grid of 100x100x100 you will have a million equation that you need to solve and that is very difficult. This is a very common problem faced in fluid dynamics. You should get in touch with one of the experts in these field.