一阶微分方程的常数变易法/洛谷P6613

一阶微分方程的常数变易法

(1) 一阶齐次线性微分方程

\[\begin{aligned} F'(x)&=P(x)F(x)\\ \dfrac{1}{F(x)}\times F'(x)&=P(x)\\ (\ln F(x))'&=P(x)\\ \ln F(x)&=\int P(x) \text dx+\ln C\\ F(x)&= Ce^{\int P(x) \text dx}\\ \end{aligned}\]

(2) 一阶非齐次线性微分方程

\[F'(x)=P(x)F(x)+Q(x) \]

考虑拆成两个只与 \(P(x)\) 或 \(Q(x)\) 其中一个相关的方程，分别求解后代回去计算 \(F(x)\)。

设 \(F(x)=u(x)v(x)\)

\[\begin{aligned} u'(x)v(x)+u(x)v'(x)&=P(x)u(x)v(x)+Q(x)\\ u'(x)v(x)+u(x)(v'(x)-P(x)v(x))&=Q(x)\\ \end{aligned}\]

若 \(v'(x)-P(x)v(x)=0\) 就可以消掉 \(P(x)\) 求解了。由 (1) 得此时 \(v(x)=C_ve^{\int P(x) \text dx}\)。带入原方程

\[\begin{aligned} u'(x)v(x)&=Q(x)\\ u'(x)C_ve^{\int P(x) \text dx}&=Q(x)\\ u'(x)&=\dfrac{Q(x)}{C_ve^{\int P(x) \text dx}}\\ u(x)&=C_v^{-1}(\int Q(x)e^{-\int P(x) \text dx}\text dx+C_u)\\ \end{aligned}\]

\(\displaystyle\therefore F(x)=u(x)v(x)=(\int Q(x)e^{-\int P(x) \text dx}\text dx+C_u)e^{\int P(x) \text dx}\)

发现解的形式与对应齐次方程的解 $ F(x)= Ce^{\int P(x) \text dx}$ 相似，只是把常数 \(C\) 换成了关于 \(x\) 的多项式 \(\displaystyle C(x)=\int Q(x)e^{-\int P(x) \text dx}\text dx+C_u\)

所以可以得到微分方程的常数变易法：先解对应的齐次方程，将常数换为待定多项式后代回原方程求解。

(3) 一阶非线性微分方程

\[\frac{\text dF(x)}{\text dx} = A(x)\text e^{F(x)-1}+B(x) \]

考虑对应的齐次方程

\[\begin{aligned} \dfrac{\text dF(x)}{\text dx}&=A(x)e^{F(x)-1}\\ \int\dfrac{\text dF(x)}{e^{F(x)-1}}&=\int A(x)\text dx\\ -e^{1-F(x)}&=\int A(x)\text dx-C\\ 1-F(x)&=\ln(C-\int A(x)\text dx)\\ F(x)&=1-\ln(C-\int A(x)\text dx)\\ \end{aligned}\]

设原方程的解为 \(\displaystyle F(x)=1-\ln(C(x)-\int A(x)\text dx)\)，代入得

\[\begin{aligned} (1-\ln(C(x)-\int A(x)\text dx))'&=\dfrac{A(x)}{\displaystyle C(x)-\int A(x)\text dx}+B(x)\\ -\dfrac{C'(x)-A(x)}{\displaystyle C(x)-\int A(x)\text dx}&=\dfrac{A(x)}{\displaystyle C(x)-\int A(x)\text dx}+B(x)\\ C'(x)&=-B(x)(C(x)-\int A(x)\text dx)\\ &=-B(x)C(x)+B(x)\int A(x)\text dx\\ \end{aligned}\]

由 (2)
得

\[\displaystyle C(x)=(\int (\int A(x)\text dx )B(x)e^{\int B(x) \text dx}\text dx+C_0)e^{-\int B(x) \text dx}\]

\(\therefore\)

\[\begin{aligned} F(x)&=1-\ln(C(x)-\int A(x)\text dx)\\ &=1-\ln((\int (\int A(x)\text dx)B(x)e^{\int B(x) \text dx}\text dx+C_0)e^{-\int B(x) \text dx}-\int A(x)\text dx)\\ &=1+\int B(x)\text dx-\ln(\int (\int A(x)\text dx)B(x)e^{\int B(x) \text dx}\text dx+C_0-\int A(x)\text dxe^{\int B(x) \text dx})\\ \end{aligned}\]

取 \(C_0=1\) 即可。

代码

#include<bits/stdc++.h>
#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
#define bceil(n) (1<<__lg(n-1)+1)
using namespace std;
const int siz=1<<18;char buf[siz],*p1=buf,*p2=buf,obuf[siz],*p3=obuf;
#define getchar() p1==p2&&(p2=(p1=buf)+fread(buf,1,siz,stdin),p1==p2)?EOF:*p1++
#define putchar(x) (p3-obuf<siz)?(*p3++=x):(fwrite(obuf,p3-obuf,1,stdout),p3=obuf,*p3++=x)
int read(){
    int a=0;char ch=getchar();
    while(ch<'0'||ch>'9') ch=getchar();
    while(ch>='0'&&ch<='9') a=(a<<3)+(a<<1)+(ch^'0'),ch=getchar();
    return a;
} 
void write(int a){
    if(a>9) write(a/10); 
    putchar(a%10+'0');
}
const int MAXN=1e6+10,P=998244353,G=3,Gi=332748118,Img=86583718;
int l,r[MAXN],inv[MAXN],sav[MAXN<<1];
ll qpow(ll a,ll b=P-2){
    if(a==1) return 1;
    ll ans=1;
    while(b){if(b&1) ans=ans*a%P;a=a*a%P;b>>=1;}
    return ans;
}
void tpre(int lim){
    if(l==lim) return;l=lim;
    for(int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)?lim>>1:0);
}
void px(int *A,int *B,int n){for(int i=0;i<n;i++) A[i]=1ll*A[i]*B[i]%P;} 
void NTT(int *A,int lim,int type){
    tpre(lim);
    static ull f[MAXN<<1],w[MAXN];w[0]=1;
    for(int i=0;i<lim;i++) f[i]=(((ll)P<<5)+A[r[i]])%P;
    for(int mid=1;mid<lim;mid<<=1){
        ull Wn=qpow(type+1?G:Gi,(P-1)/(mid+mid));
        for(int i=1;i<mid;i++)w[i]=w[i-1]*Wn%P;
        for(int j=0;j<lim;j+=mid+mid){
            for(int k=0;k<mid;k++){
                int x=w[k]*f[j|mid|k]%P;
                f[j|mid|k]=f[j|k]+P-x;
                f[j|k]+=x;
            }   
        }if(mid==(1<<10)){for(int i=0;i<lim;i++) f[i]%=P;}
    }if(type-1){
        ull inv=qpow(lim);
        for(int i=0;i<lim;i++) A[i]=f[i]%P*inv%P;
    }else for(int i=0;i<lim;i++) A[i]=f[i]%P;
}
void mul(int *A,int *B,int la,int lb){
    int lim=bceil(la+la);
    cpy(sav,B,lim);clr(sav+la,lim-la);
    NTT(A,lim,1);NTT(sav,lim,1);
    px(A,sav,lim);NTT(A,lim,-1);
    clr(A+lb,lim-lb);clr(sav,lim);
} 
void invp(int *A,int lim){
    int n=bceil(lim);
    static int w[MAXN<<1],r[MAXN<<1];
    w[0]=qpow(A[0]);
    for (int ln=2;ln<=n;ln<<=1){
        for(int i=0;i<(ln>>1);i++) r[i]=w[i];
        cpy(sav,A,ln);NTT(sav,ln,1);NTT(r,ln,1);px(r,sav,ln);
        NTT(r,ln,-1);clr(r,ln>>1);cpy(sav,w,ln);NTT(sav,ln,1);
        NTT(r,ln,1);px(r,sav,ln);NTT(r,ln,-1);
        for(int i=ln>>1;i<ln;i++) w[i]=(w[i]*2ll-r[i]+P)%P;
    }cpy(A,w,lim);clr(sav,n);clr(w,n);clr(r,n);
}
void deriv(int *A,int lim){
    for(int i=1;i<lim;i++) A[i-1]=1ll*A[i]*i%P;
    A[lim-1]=0;
}
void inv_init(int lim){
    inv[1]=1;
    for(int i=2;i<=lim;i++) inv[i]=1ll*inv[P%i]*(P-P/i)%P;
}
void integ(int *A,int lim){
    for(int i=lim;i;i--) A[i]=1ll*A[i-1]*inv[i]%P;
    A[0]=0;
}
void lnp(int *A,int lim){
    static int w[MAXN<<1];
    cpy(w,A,lim);
    invp(w,lim);deriv(A,lim);
    mul(A,w,lim,lim);
    integ(A,lim-1);
    clr(w,lim);
}
void exp(int *A,int lim){
    static int s[MAXN<<1],s2[MAXN<<1];
    int n=bceil(lim);
    s2[0]=1;
    for(int ln=2;ln<=n;ln<<=1){
        cpy(s,s2,ln>>1);lnp(s,ln);
        for(int i=0;i<ln;i++) s[i]=(A[i]-s[i]+P)%P;
        s[0]=(s[0]+1)%P;
        mul(s2,s,ln,ln);
    }cpy(A,s2,lim);clr(s,n);clr(s2,n);
}
int n,m,A[MAXN],B[MAXN],F[MAXN],s[MAXN];
int main(){
	n=read()+1;inv_init(n);
	for(int i=0;i<n;i++) A[i]=read();
	for(int i=0;i<n;i++) B[i]=read();
	integ(A,n);cpy(F,B,n);integ(F,n);
	cpy(s,F,n);exp(s,n);mul(s,A,n,n);
	mul(B,s,n,n);integ(B,n);B[0]++;
	for(int i=0;i<n;i++) B[i]=(B[i]-s[i]+P)%P;
	lnp(B,n);F[0]++;
	for(int i=0;i<n;i++) F[i]=(F[i]-B[i]+P)%P;
	for(int i=0;i<n;i++) write(F[i]),putchar(' ');
    fwrite(obuf,p3-obuf,1,stdout);
    return 0;
}