C++代码实现MATLAB中的pade函数功能

// pade.cpp// 用 C++ 实现 MATLAB 的 pade 函数（幂级数版本 + 纯延时版本）// 编译： g++ -std=c++17 -O2 pade.cpp -o pade#include<vector>#include<complex>#include<iostream>#include<iomanip>#include<cmath>#include<stdexcept>#include<algorithm>#include<string>namespacempad{usingComplex=std::complex<double>;usingVecC=std::vector<Complex>;usingMatC=std::vector<VecC>;// ------------------------------------------------------------------// 列主元高斯消元求解 A*x = b// ------------------------------------------------------------------VecCsolveLinearSystem(constMatC&A_in,constVecC&b_in){constintn=static_cast<int>(A_in.size());if(n==0)return{};MatC A=A_in;VecC b=b_in;for(inti=0;i<n;++i)A[i].push_back(b[i]);// 增广矩阵for(intc=0;c<n;++c){intpiv=c;doublebest=std::abs(A[c][c]);for(intr=c+1;r<n;++r){if(std::abs(A[r][c])>best){best=std::abs(A[r][c]);piv=r;}}if(best<1e-14)throwstd::runtime_error("pade: 线性方程组奇异，无法求得 Pade 系数");if(piv!=c)std::swap(A[c],A[piv]);for(intr=c+1;r<n;++r){Complex f=A[r][c]/A[c][c];for(intk=c;k<=n;++k)A[r][k]-=f*A[c][k];}}VecCx(n);for(intr=n-1;r>=0;--r){Complex s=A[r][n];for(intc=r+1;c<n;++c)s-=A[r][c]*x[c];x[r]=s/A[r][r];}returnx;}// ------------------------------------------------------------------// 核心：由幂级数系数（升幂）c[0..L+M] 求 [L/M] 阶 Pade 逼近// 输出 P（分子系数，升幂，长度 L+1）// Q（分母系数，升幂，长度 M+1，Q[0] == 1）// ------------------------------------------------------------------voidpadeLM(constVecC&c,intL,intM,VecC&P,VecC&Q){constintN=L+M;if(static_cast<int>(c.size())<N+1)throwstd::runtime_error("pade: 幂级数系数个数不足 (需要 >= L+M+1)");// 解 Q(x) = 1 + q1 x + ... + qM x^M// 匹配方程 (k = L+1 .. L+M): 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(intr=0;r<M;++r){constintk=L+1+r;for(intj=0;j<M;++j){constintidx=k-(j+1);A[r][j]=(idx>=0)?c[idx]:Complex(0.0,0.0);}rhs[r]=-c[k];}Q.assign(M+1,Complex(0.0,0.0));Q[0]=Complex(1.0,0.0);if(M>0){VecC q=solveLinearSystem(A,rhs);for(intj=0;j<M;++j)Q[j+1]=q[j];}// 回代 P(x) = p0 + p1 x + ... + pL x^L// p_k = sum_{i=0..min(k,M)} q_i * c_{k-i}P.assign(L+1,Complex(0.0,0.0));for(intk=0;k<=L;++k){Complexs(0.0,0.0);constintimax=std::min(k,M);for(inti=0;i<=imax;++i)s+=Q[i]*c[k-i];P[k]=s;}}// 升幂 -> 降幂（MATLAB 传递函数约定）std::vector<double>toDescending(constVecC&a){constintn=static_cast<int>(a.size());std::vector<double>out(n);for(inti=0;i<n;++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::vector<double>&coeffsAscending,intN,std::vector<double>&num,std::vector<double>&den){if(N<0)throwstd::runtime_error("pade: N 必须为非负整数");constintL=N/2;constintM=N-L;VecCc(coeffsAscending.begin(),coeffsAscending.end());VecC P,Q;padeLM(c,L,M,P,Q);num=toDescending(P);den=toDescending(Q);}// ------------------------------------------------------------------// [num,den] = pade(T, N) —— exp(-T*s) 的 N 阶 Pade 逼近// ------------------------------------------------------------------voidpadeDelay(doubleT,intN,std::vector<double>&num,std::vector<double>&den){if(N<0)throwstd::runtime_error("pade: N 必须为非负整数");constintL=N/2;constintM=N-L;// exp(-T*s) 在 s=0 的泰勒系数 c_k = (-T)^k / k!std::vector<double>c(L+M+1);c[0]=1.0;for(intk=1;k<=L+M;++k)c[k]=c[k-1]*(-T)/static_cast<double>(k);pade(c,N,num,den);}voidprintVector(conststd::string&name,conststd::vector<double>&v){std::cout<<name<<" = [";for(size_t i=0;i<v.size();++i){std::cout<<v[i];if(i+1<v.size())std::cout<<", ";}std::cout<<"]\n";}}// namespace mpad// ==================================================================intmain(){usingnamespacempad;std::cout<<std::setprecision(10);// 例 1：exp(x) 的 [2/3] Pade// coeffs = 1./factorial(0:5)std::vector<double>coeffs={1.0,1.0,0.5,1.0/6.0,1.0/24.0,1.0/120.0};std::vector<double>num,den;pade(coeffs,5,num,den);std::cout<<"--- [2/3] Pade of exp(x) ---\n";printVector("num",num);printVector("den",den);// 例 2：exp(-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);// 例 3：exp(-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;}