MacCormack差分格式的全局误差分析
MacCormack差分格式的全局誤差分析
問題
請分析MacCormack差分格式的全局誤差:
un+12=un+hun′u_{n+\frac{1}{2}} = u_n + h u'_nun+21??=un?+hun′?
un+1=12(un+un+12+hun+12′)u_{n+1} =\frac{1}{2} \left( u_n + u_{n+\frac{1}{2}} + h u'_{n+\frac{1}{2}} \right)un+1?=21?(un?+un+21??+hun+21?′?)
解答思路
首先,得把這個格式的寫法改改,寫成:
(MC)u~n+1=un+hun′un+1=12[un+u~n+1+hu~n+1′]\begin{aligned} \tilde{u}_{n+1} &=u_{n}+h u_{n}^{\prime} \\ u_{n+1} &=\frac{1}{2}\left[u_{n}+\tilde{u}_{n+1}+h \tilde{u}_{n+1}^{\prime}\right] \end{aligned} \tag{MC} u~n+1?un+1??=un?+hun′?=21?[un?+u~n+1?+hu~n+1′?]?(MC)
或者
u~n+12=un+hun′\tilde u_{n+\frac{1}{2}} = u_n + h u'_nu~n+21??=un?+hun′?
un+1=12(un+u~n+12+hu~n+12′)u_{n+1} =\frac{1}{2} \left( u_n +\tilde u_{n+\frac{1}{2}} + h \tilde u'_{n+\frac{1}{2}} \right)un+1?=21?(un?+u~n+21??+hu~n+21?′?)
因為涉及到位移算子EEE移一格還是移半格的問題,不然會出問題的,我覺得是這樣,但也可能不對。下面就以改成第一種方式來算。
對于代表性O(shè)DE
dudt=u′=λu+aeμt\frac{d u}{d t}=u^{\prime}=\lambda u+a e^{\mu t} dtdu?=u′=λu+aeμt
它的精確解是:
u(nh)=c(eλh)n+a(e(μh))nμ?λu(n h)=c\left(e^{\lambda h}\right)^{n}+\frac{a\left(e^{(\mu h)}\right)^{n}}{\mu-\lambda} u(nh)=c(eλh)n+μ?λa(e(μh))n?
借助位移算子,將ODE代入到數(shù)值格式MCMCMC,
整理理一下可以得到:
P(E)un=Q(E)?aeμhnP(E) u_{n}=Q(E) \cdot a e^{\mu h n} P(E)un?=Q(E)?aeμhn
那么數(shù)值解實際上就是:
un=∑k=1Kck(σk)n+aeμhn?Q(eμh)P(eμh)u_{n}=\sum_{k=1}^{K} c_{k}\left(\sigma_{k}\right)^{n}+a e^{\mu h n} \cdot \frac{Q\left(e^{\mu h}\right)}{P\left(e^{\mu h}\right)} un?=k=1∑K?ck?(σk?)n+aeμhn?P(eμh)Q(eμh)?
我們可以通過比較數(shù)值解和精確解(泰勒展開)來得到誤差。誤差項分為特解項(第一部分)的誤差,和非特解項(第二部分)的誤差。如果數(shù)值解第一部分的σ\sigmaσ是復數(shù)的話,要分解計算模的誤差和復角的誤差。
具體計算
MCMCMC格式對于代表性O(shè)DE
dudt=u′=λu+aeμt\frac{d u}{d t}=u^{\prime}=\lambda u+a e^{\mu t} dtdu?=u′=λu+aeμt
的離散格式為:
u~n+1?(1+λh)un=aheμhn?12(1+λh)u~n+1+un+1?12un=12aheμh(n+1)\begin{aligned} \tilde{u}_{n+1}-(1+\lambda h) u_{n} &=a h e^{\mu h n} \\-\frac{1}{2}(1+\lambda h) \tilde{u}_{n+1}+u_{n+1}-\frac{1}{2} u_{n} &=\frac{1}{2} a h e^{\mu h(n+1)} \end{aligned} u~n+1??(1+λh)un??21?(1+λh)u~n+1?+un+1??21?un??=aheμhn=21?aheμh(n+1)?
寫成位移算子的形式:
Eu~n?(1+λh)un=aheμhn?12(1+λh)Eu~n+Eun?12un=12Eaheμh(n)\begin{aligned} E \tilde{u}_{n}-(1+\lambda h) u_{n} &=a h e^{\mu h n} \\-\frac{1}{2}(1+\lambda h) E \tilde{u}_{n}+E u_{n}-\frac{1}{2} u_{n} &=\frac{1}{2} E a h e^{\mu h(n)} \end{aligned} Eu~n??(1+λh)un??21?(1+λh)Eu~n?+Eun??21?un??=aheμhn=21?Eaheμh(n)?
整理成矩陣表達為:
[E?(1+λh)?12(1+λh)EE?12][u~u]n=h[112E]aeμhn(19)\left[ \begin{array}{cc}{E} & {-(1+\lambda h)} \\ {-\frac{1}{2}(1+\lambda h) E} & {E-\frac{1}{2}}\end{array}\right] \left[ \begin{array}{c}{\tilde{u}} \\ {u}\end{array}\right]_{n}=h \left[ \begin{array}{c}{1} \\ {\frac{1}{2} E}\end{array}\right] a e^{\mu h n}(19) [E?21?(1+λh)E??(1+λh)E?21??][u~u?]n?=h[121?E?]aeμhn(19)
根據(jù)克拉默法則,求P(E)P(E)P(E)和Q(E)Q(E)Q(E),得到:
P(E)=det?[E?(1+λh)?12(1+λh)EE?12]=E(E?1?λh?12λ2h2)Q(E)=det?[Eh?12(1+λh)E12hE]=12hE(E+1+λh)\begin{aligned} P(E) &=\operatorname{det} \left[ \begin{array}{cc}{E} & {-(1+\lambda h)} \\ {-\frac{1}{2}(1+\lambda h) E} & {E-\frac{1}{2}}\end{array}\right] \\=& E\left(E-1-\lambda h-\frac{1}{2} \lambda^{2} h^{2}\right) \\ Q(E) &=\operatorname{det} \left[ \begin{array}{cc}{E} & {h} \\ {-\frac{1}{2}(1+\lambda h) E} & {\frac{1}{2} h E}\end{array}\right] \\ &=\frac{1}{2} h E(E+1+\lambda h) \end{aligned} P(E)=Q(E)?=det[E?21?(1+λh)E??(1+λh)E?21??]E(E?1?λh?21?λ2h2)=det[E?21?(1+λh)E?h21?hE?]=21?hE(E+1+λh)?
特征多項式P(σ)P(\sigma)P(σ)為:
P(σ)=σ(σ?1?λh?12λ2h2)=0P(\sigma)=\sigma\left(\sigma-1-\lambda h-\frac{1}{2} \lambda^{2} h^{2}\right)=0 P(σ)=σ(σ?1?λh?21?λ2h2)=0
將其根以及PPP和QQQ代入到數(shù)值解的表達式:
un=∑k=1Kck(σk)n+aeμhn?Q(eμh)P(eμh)u_{n}=\sum_{k=1}^{K} c_{k}\left(\sigma_{k}\right)^{n}+a e^{\mu h n} \cdot \frac{Q\left(e^{\mu h}\right)}{P\left(e^{\mu h}\right)} un?=k=1∑K?ck?(σk?)n+aeμhn?P(eμh)Q(eμh)?
得到數(shù)值解的一個表達:
un=c1(1+λh+12λ2h2)n+aeμhn?12h(eμh+1+λh)eμh?1?λh?12λ2h2\begin{array}{c}{u_{n}=c_{1}\left(1+\lambda h+\frac{1}{2} \lambda^{2} h^{2}\right)^{n}+} \\ {a e^{\mu h n} \cdot \frac{\frac{1}{2} h\left(e^{\mu h}+1+\lambda h\right)}{e^{\mu h}-1-\lambda h-\frac{1}{2} \lambda^{2} h^{2}}}\end{array} un?=c1?(1+λh+21?λ2h2)n+aeμhn?eμh?1?λh?21?λ2h221?h(eμh+1+λh)??
考慮這個代表問題的精確解:
u(nh)=c(eλh)n+aeμnhμ?λu(n h)=c\left(e^{\lambda h}\right)^{n}+\frac{ae^{\mu nh}}{\mu-\lambda} u(nh)=c(eλh)n+μ?λaeμnh?
下面來考慮全局誤差,λ,h\lambda,hλ,h都為實數(shù),假定最后的時間為TTT,時間步為NNN,即T=NhT=NhT=Nh。
那么全局transient的誤差可以表示為:
Erλ≡eλT?(σ1(λh))N=eλNh?(1+λh+12λ2h2)NE r_{\lambda} \equiv e^{\lambda T}-\left(\sigma_{1}(\lambda h)\right)^{N}=e^{\lambda Nh}-\left(1+\lambda h+\frac{1}{2} \lambda^{2} h^{2}\right)^{N} Erλ?≡eλT?(σ1?(λh))N=eλNh?(1+λh+21?λ2h2)N
特解的全局誤差定義為第二部的比值和1的差值:
Eμ≡(μ?λ)Q(eμh)P(eμh)?1=(μ?λ)12h(eμh+1+λh)eμh?1?λh?12λ2h2?1E_{\boldsymbol{\mu}} \equiv(\mu-\lambda) \frac{Q\left(e^{\mu h}\right)}{P\left(e^{\mu h}\right)}-1={(\mu - \lambda)\frac{\frac{1}{2} h\left(e^{\mu h}+1+\lambda h\right)}{e^{\mu h}-1-\lambda h-\frac{1}{2} \lambda^{2} h^{2}}}-1 Eμ?≡(μ?λ)P(eμh)Q(eμh)??1=(μ?λ)eμh?1?λh?21?λ2h221?h(eμh+1+λh)??1
PS: 預測矯正格式,用兩個公式做誤差分析和通過替換,消掉中間變量,把兩個公式變成一個再做誤差分析,結(jié)果不盡相同。
總結(jié)
以上是生活随笔為你收集整理的MacCormack差分格式的全局误差分析的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: .NET c#中调用地图
- 下一篇: 9.6 模拟试题