一阶微分方程的常数变易法
(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)
得
\(\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;
}