Stable method for numerically solving matrix ODE I’m…
0 reactions
Stable method for numerically solving matrix ODE
I’m setting up a simulation (RCWA in electromagnetism) which requires me to solve d/dz y=Ay. However, A is a massive matrix with a large L1 norm. This makes diaganolization impractical (besides for a very crude simulation), and taking exp(A) seems not to work well (I am assuming there is floating point error with my tiny scale factor that causes exp(A/N)^N to lose a lot of accuracy). Even if I implemented some super stable algorithm I’m pretty sure I’d eventually surpass the floating point maximum making this pointless.
I will note that there is reason to believe the equation should still be solvable even with these issues—y should be a relatively nice vector, maybe with elements that are close to 0. I don’t think it’ll be close to machine epsilon though.
So now I’m 0/2 for the most common methods to solve such an equation. I am wondering if there is any other approach worth trying. I’m wondering if maybe some high order implicit ODE solver would work well. I’d also guess there may be some Krylov method for computing exp(A)x but I haven’t seen any (and would kind of prefer something that is widely implemented or won’t take a super long time to implement). I was also thinking Galerkin methods may be applicable but this seems like it may require a very fine discretization. I’d appreciate any suggestions as I’m a bit stuck.
It might be worth mentioning A (should) have a pretty decent preconditioner if this may make some options viable. Also, A is a block matrix of the form [0,P;Q,0], but P and Q don’t have a great structure (essentially Toeplitz matricies sandwiched between diaganol matricies). Otherwise there’s not much else to the problem.
#science
source
I’m setting up a simulation (RCWA in electromagnetism) which requires me to solve d/dz y=Ay. However, A is a massive matrix with a large L1 norm. This makes diaganolization impractical (besides for a very crude simulation), and taking exp(A) seems not to work well (I am assuming there is floating point error with my tiny scale factor that causes exp(A/N)^N to lose a lot of accuracy). Even if I implemented some super stable algorithm I’m pretty sure I’d eventually surpass the floating point maximum making this pointless.
I will note that there is reason to believe the equation should still be solvable even with these issues—y should be a relatively nice vector, maybe with elements that are close to 0. I don’t think it’ll be close to machine epsilon though.
So now I’m 0/2 for the most common methods to solve such an equation. I am wondering if there is any other approach worth trying. I’m wondering if maybe some high order implicit ODE solver would work well. I’d also guess there may be some Krylov method for computing exp(A)x but I haven’t seen any (and would kind of prefer something that is widely implemented or won’t take a super long time to implement). I was also thinking Galerkin methods may be applicable but this seems like it may require a very fine discretization. I’d appreciate any suggestions as I’m a bit stuck.
It might be worth mentioning A (should) have a pretty decent preconditioner if this may make some options viable. Also, A is a block matrix of the form [0,P;Q,0], but P and Q don’t have a great structure (essentially Toeplitz matricies sandwiched between diaganol matricies). Otherwise there’s not much else to the problem.
#science
source
| 0 mlx | |
| 0 | |
| 0 |