Utilisation optimale du fractionnement de Strang (pour l'équation de diffusion de la réaction)

9

J'ai fait une observation étrange en calculant la solution d'une simple équation de diffusion de réaction 1D:

ta=2x2aab

tb=ab

tc=a

La valeur initiale de est une constante ( b ( 0 , x ) = b 0 ), et je ne suis intéressé que par l'intégrale sur a de 0 à 1 ( 1 0 a ( t , x ) d t ). Le but de c et l'équation bb(0,x)=b0a0101a(t,x)dtcest juste pour évaluer cette intégrale.tc=a

J'ai utilisé un schéma de fractionnement de Strang pour le couplage entre la diffusion et la réaction (une réaction en demi-étape, puis une diffusion en étape complète, puis à nouveau une réaction en demi-étape), un schéma de Crank Nicholson pour la diffusion et une solution analytique pour la réaction ( y compris l'équation ).tc=a

Parce qu'une étape de la solution analytique était plus d'un facteur 3 plus lente qu'une étape du schéma de Crank Nicholson, j'ai essayé de faire plus d'une étape de Crank Nicholson pour chaque étape de réaction. J'espérais m'en sortir avec moins d'étape du schéma de fractionnement Strang, afin d'être plus rapide dans l'ensemble.

Cependant, l'effet inverse peut être observé, à savoir que beaucoup plus d'étapes pour le schéma de fractionnement de Strang sont nécessaires si plus d'une étape Crank Nicholson est utilisée. (Je ne m'intéresse qu'à la précision de l'intégrale sur , qui semble converger plus vite que a lui - même.) Après m'être demandé pendant un certain temps, j'ai remarqué que le même effet se produit également pour b ( t , x ) = b 0 = 0 , et que je comprends même pourquoi pour ce cas. Le fait est que si je fais exactement un pas de Crank Nicholson, alors le schéma global se transforme en une règle trapézoïdale (si b ( t ) = 0 ).aab(t,x)=b0=0b(t)=0

Donc si je traitais dans le cadre de l'étape de diffusion, l'augmentation du nombre d'étapes Crank Nicholson (probablement) ne conduirait pas à une précision globale réduite (comme observé). Mais cela semble aller à l'encontre de l'objectif d'utiliser une solution analytique pour la partie réactionnelle (non linéaire et potentiellement très rigide) du système.tc=a

Voici donc ma question: existe-t-il une meilleure façon de traiter dans le cadre d'un fractionnement de Strang, que de le traiter dans le cadre de l'étape de réaction, ou de le traiter dans le cadre de l'étape de diffusion. Je veux éviter d'être «forcé» d'utiliser exactement une étape Crank Nicholson pour la diffusion. (Par exemple, en 3D, je préférerais résoudre la diffusion analytiquement par une FFT au lieu d'utiliser Crank Nicholson. Bien sûr, je peux également combiner FFT avec Crank Nicholson, donc ce n'est pas si grave.)tc=a

Thomas Klimpel
la source
Dans people.maths.ox.ac.uk/dellar/OperatorLB.html , un effet similaire semble être décrit. La conclusion est qu'il est crucial d'utiliser Crank Nicholson au lieu de la solution exacte. Alors peut-être que la réponse à ma question est un simple non.
Thomas Klimpel
Quelque chose ne va pas avec vos équations. n'apparaît pas dans les deux premiers, ce qui rend le couplage unidirectionnel et signifie que vous pouvez calculer c à tout t comme étape de post-traitement. cct
Bill Barth
@BillBarth J'ai changé la question pour clarifier le rôle de . Donc c est juste un moyen de calculer 1 0 a ( t , x ) d t . Veuillez me faire savoir si vous avez des suggestions sur la façon de calculer cette intégrale plus précise (que ce que j'obtiens de la combinaison de la séparation de Strang et de Crank Nicholson décrite ci-dessus), en utilisant potentiellement une étape de post-traitement. cc01a(t,x)dt
Thomas Klimpel
Cela fait maintenant longtemps, mais avez-vous reconnu que ce système d'équations peut être écrit comme une PDE parabolique en avec un terme de réaction exponentielle? Je suppose que je me demande si vous voulez vraiment résoudre ce système à 3 variables au lieu d'un système simplifié. c
Bill Barth
@BillBarth Je serais intéressé d'apprendre comment ce système peut être écrit comme un PDE parabolique avec un terme de réaction exponentielle. La vitesse de la solution de ce modèle est un facteur limitant lors de l'étalonnage du modèle (qui peut prendre plusieurs heures), même si la précision utilisée en matière d'intégration temporelle est assez loin de la pleine convergence.
Thomas Klimpel

Réponses:

6

Je vais écrire ceci comme une réponse bien qu'il ne réponde pas directement à la question.

En branchant la deuxième équation et la troisième équation dans la première, et en branchant la troisième dans la seconde, donnez ensemble: Réorganiser ces deux donne:

2ct2=2x2ct+btbt=(ct)b
Maintenant, nous pouvons intégrer ces deux éléments une fois ent, en laissant la première équation: c
t(ct2cx2b)=01b(bt)=ct
t
ct2cx2=b+A(x)
A(x)=a02c0x2b0 Cela conduit à la solution lnb(x,t)-lnb0(x)=-c(x)+c0(x) Ou lnb
0t1b(x,t)(b(x,t)t)dt=0tc(x,t)tdt
lnb(x,t)lnb0(x)=c(x)+c0(x)
Exponentiating donne: b=b0ec0-c Et, enfin, le brancher dans le PDE pourcdonne c
lnbb0=c+c0
b=b0ec0c
c
ct2cx2=b0ec0c+A(x)

ccc0c0(x)=0

ct=2cx2+a0(1ec)b0
c>c0

b0=0

t=1c(x,1)

Bill Barth
la source
ca
Aucun problème. Je me moquais de moi après notre premier échange de commentaires. J'espère que c'est utile.
Bill Barth