// pade.cpp// 用 C 实现 MATLAB 的 pade 函数幂级数版本 纯延时版本// 编译 g -stdc17 -O2 pade.cpp -o pade#includevector#includecomplex#includeiostream#includeiomanip#includecmath#includestdexcept#includealgorithm#includestringnamespacempad{usingComplexstd::complexdouble;usingVecCstd::vectorComplex;usingMatCstd::vectorVecC;// ------------------------------------------------------------------// 列主元高斯消元求解 A*x b// ------------------------------------------------------------------VecCsolveLinearSystem(constMatCA_in,constVecCb_in){constintnstatic_castint(A_in.size());if(n0)return{};MatC AA_in;VecC bb_in;for(inti0;in;i)A[i].push_back(b[i]);// 增广矩阵for(intc0;cn;c){intpivc;doublebeststd::abs(A[c][c]);for(intrc1;rn;r){if(std::abs(A[r][c])best){beststd::abs(A[r][c]);pivr;}}if(best1e-14)throwstd::runtime_error(pade: 线性方程组奇异无法求得 Pade 系数);if(piv!c)std::swap(A[c],A[piv]);for(intrc1;rn;r){Complex fA[r][c]/A[c][c];for(intkc;kn;k)A[r][k]-f*A[c][k];}}VecCx(n);for(intrn-1;r0;--r){Complex sA[r][n];for(intcr1;cn;c)s-A[r][c]*x[c];x[r]s/A[r][r];}returnx;}// ------------------------------------------------------------------// 核心由幂级数系数升幂c[0..LM] 求 [L/M] 阶 Pade 逼近// 输出 P分子系数升幂长度 L1// Q分母系数升幂长度 M1Q[0] 1// ------------------------------------------------------------------voidpadeLM(constVecCc,intL,intM,VecCP,VecCQ){constintNLM;if(static_castint(c.size())N1)throwstd::runtime_error(pade: 幂级数系数个数不足 (需要 LM1));// 解 Q(x) 1 q1 x ... qM x^M// 匹配方程 (k L1 .. LM): q1*c[k-1] ... qM*c[k-M] -c[k]MatCA(M,VecC(M,Complex(0.0,0.0)));VecCrhs(M,Complex(0.0,0.0));for(intr0;rM;r){constintkL1r;for(intj0;jM;j){constintidxk-(j1);A[r][j](idx0)?c[idx]:Complex(0.0,0.0);}rhs[r]-c[k];}Q.assign(M1,Complex(0.0,0.0));Q[0]Complex(1.0,0.0);if(M0){VecC qsolveLinearSystem(A,rhs);for(intj0;jM;j)Q[j1]q[j];}// 回代 P(x) p0 p1 x ... pL x^L// p_k sum_{i0..min(k,M)} q_i * c_{k-i}P.assign(L1,Complex(0.0,0.0));for(intk0;kL;k){Complexs(0.0,0.0);constintimaxstd::min(k,M);for(inti0;iimax;i)sQ[i]*c[k-i];P[k]s;}}// 升幂 - 降幂MATLAB 传递函数约定std::vectordoubletoDescending(constVecCa){constintnstatic_castint(a.size());std::vectordoubleout(n);for(inti0;in;i)out[i]a[n-1-i].real();returnout;}// ------------------------------------------------------------------// [num,den] pade(coeffs, N)// coeffs : f(x) 幂级数系数升幂// N : 总阶数L floor(N/2)M N - L// ------------------------------------------------------------------voidpade(conststd::vectordoublecoeffsAscending,intN,std::vectordoublenum,std::vectordoubleden){if(N0)throwstd::runtime_error(pade: N 必须为非负整数);constintLN/2;constintMN-L;VecCc(coeffsAscending.begin(),coeffsAscending.end());VecC P,Q;padeLM(c,L,M,P,Q);numtoDescending(P);dentoDescending(Q);}// ------------------------------------------------------------------// [num,den] pade(T, N) —— exp(-T*s) 的 N 阶 Pade 逼近// ------------------------------------------------------------------voidpadeDelay(doubleT,intN,std::vectordoublenum,std::vectordoubleden){if(N0)throwstd::runtime_error(pade: N 必须为非负整数);constintLN/2;constintMN-L;// exp(-T*s) 在 s0 的泰勒系数 c_k (-T)^k / k!std::vectordoublec(LM1);c[0]1.0;for(intk1;kLM;k)c[k]c[k-1]*(-T)/static_castdouble(k);pade(c,N,num,den);}voidprintVector(conststd::stringname,conststd::vectordoublev){std::coutname [;for(size_t i0;iv.size();i){std::coutv[i];if(i1v.size())std::cout, ;}std::cout]\n;}}// namespace mpad// intmain(){usingnamespacempad;std::coutstd::setprecision(10);// 例 1exp(x) 的 [2/3] Pade// coeffs 1./factorial(0:5)std::vectordoublecoeffs{1.0,1.0,0.5,1.0/6.0,1.0/24.0,1.0/120.0};std::vectordoublenum,den;pade(coeffs,5,num,den);std::cout--- [2/3] Pade of exp(x) ---\n;printVector(num,num);printVector(den,den);// 例 2exp(-s) 的 [2/2] Pade等价 pade(1,4)padeDelay(1.0,4,num,den);std::cout\n--- [2/2] Pade of exp(-s) ---\n;printVector(num,num);printVector(den,den);// 例 3exp(-0.1s) 的 [1/2] Pade等价 pade(0.1,3)padeDelay(0.1,3,num,den);std::cout\n--- [1/2] Pade of exp(-0.1s) ---\n;printVector(num,num);printVector(den,den);return0;}