一阶微分方程

· · 算法·理论

常数变易法

考虑以下一阶线性微分方程

\dfrac{\text{d}y}{\text{d}x}=P(x)y+Q(x)

先解齐次方程

\begin{aligned} \dfrac{\text{d}y}{\text{d}x}&=P(x)y\\ \dfrac{\text{d}y}{y}&=P(x)\text{d}x\\ \ln y&=P_1(x)+C\\ y&=C\exp(P_1(x)) \end{aligned}

其中

P_1(x)=\int P(x)\text{d}x

根据常数变易法，对于原非齐次方程有

y=C(x)\exp(P_1(x))

即将 C 改为关于 x 的函数 C(x)。

将 y 代回原方程即可解出 C(x)。

这个方法可以拓展到非线性微分方程。

P6613

考虑以下方程

\frac{\text dF(x)}{\text dx} \equiv A(x)\text e^{F(x)-1}+B(x) \pmod{x^n}

设 y=F(x)，先解齐次方程

\begin{aligned} \dfrac{\text dy}{\text dx}&=A(x)\exp(y-1)\\ \dfrac{\text{d}y}{\exp(y-1)}&=A(x)\text{d}x\\ \exp(1-y)\text{d}y&=A(x)\text{d}x\\ -\exp(1-y)&=A_1(x)-C\\ \exp(1-y)&=C-A_1(x)\\ 1-y&=\ln(C-A_1(x))\\ y&=1-\ln(C-A_1(x))\\ y&=1+\ln\dfrac{1}{C-A_1(x)} \end{aligned}

其中

A_1(x)=\int A(x)\text{d}x

根据常数变易法，对于原非齐次方程有

\begin{aligned} y&=1-\ln(C(x)-A_1(x))\\ y&=1+\ln\dfrac{1}{C(x)-A_1(x)} \end{aligned}

代回原方程

\begin{aligned} \dfrac{\text dy}{\text dx}&=A(x)\exp(y-1)+B(x)\\ -\dfrac{C'(x)-A(x)}{C(x)-A_1(x)}&=\dfrac{A(x)}{C(x)-A_1(x)}+B(x)\\ \dfrac{C'(x)}{C(x)-A_1(x)}&=-B(x)\\ C'(x)&=-B(x)C(x)+A_1(x)B(x) \end{aligned}

设 z=C(x)

z'=-B(x)z+A_1(x)B(x)

解此方程，求出通解

z=C_1(x)\exp(-B_1(x))

其中

B_1(x)=\int B(x)\text{d}x

代入得

\begin{aligned} z'&=-B(x)z+A_1(x)B(x)\\ C_1'(x)\exp(-B_1(x))-C_1(x)\exp(-B_1(x))B(x)&=-B(x)C_1(x)\exp(-B_1(x))+A_1(x)B(x)\\ C_1'(x)\exp(-B_1(x))&=A_1(x)B(x)\\ C_1(x)&=\int A_1(x)B(x)\exp(B_1(x))\text{d}x+C \end{aligned}

继续回代

C(x)=\left(\int A_1(x)B(x)\exp(B_1(x))\text{d}x+C\right)\exp(-B_1(x)) \begin{aligned} F(x)&=1-\ln\left(\left(\int A_1(x)B(x)\exp(B_1(x))\text{d}x+C\right)\exp(-B_1(x))-A_1(x)\right)\\ &=1+B_1(x)-\ln\left(\int A_1(x)B(x)\exp(B_1(x))\text{d}x-A_1(x)\exp(B_1(x))+C\right) \end{aligned}

当 x=0，F(x)=1，故取 C=1。

注意到

B(x)\exp(B_1(x))=\exp(B_1(x))'

以及

A_1(x)\exp(B_1(x))'=(A_1(x)\exp(B_1(x)))'-A(x)\exp(B_1(x))

所以

F(x)=1+B_1(x)-\ln(1-\int {A(x)\exp(B_1(x))\text{d}x})

直接计算即可。

Code

int main(){
    int n;
    scanf("%d",&n);
    n++;
    for(int i=0;i<=n-1;i++){
        scanf("%llu",&a[i]);
    }
    for(int i=0;i<=n-1;i++){
        scanf("%llu",&b[i]);
        b1[i]=b[i];
    }
    inv[1]=1;
    for(int i=2;i<=n;i++){
        inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    }
    ITG(b1,n);
    EXP(b1,c,n);
    CPY(ans,b1,n);
    ans[0]=(ans[0]+1)%mod;
    CLR(b1,n);
    MUL(a,c,n<<1);
    CLR(a+n,n);
    ITG(a,n);
    a[0]=(1-a[0]+mod)%mod;
    for(int i=1;i<=n-1;i++){
        a[i]=(mod-a[i])%mod;
    }
    LN(a,c,n);
    DEC(ans,c,n);
    for(int i=0;i<=n-1;i++){
        printf("%llu ",ans[i]);
    }
}