如何将 Prony 系列与 PGFPlot 中的数据进行拟合？

Question

我正在尝试将 Prony 序列与 PGFPlots 中的粘弹性数据进行拟合。我遇到的问题是如何以数学方式定义能够gnuplot识别函数并进行拟合的 Prony 序列。Prony 序列是\sum_{i=0}^{k} G_k * \frac{x*a_k}{1+(x*a_k)^2}。有人知道怎么做吗？这是我使用的完整代码：

\documentclass{standalone}

\usepackage{pgfplots}
\usepackage{pgfplotstable}
\usepackage{tikz}
\usetikzlibrary{tikzmark}
\usepackage{gnuplottex}
\usepackage{filecontents}

\pgfkeys{/pgf/number format/.cd,1000 sep={}}




\pgfplotsset{compat=1.15}
\pgfplotsset{label style={font=\Large},
            tick label style={font=\Large}}

\write18{}

\begin{filecontents}{double.dat}
100     96413.1
68.1292 81638.4
46.4159 68107.2
31.6228 56025.6
21.5443 45408.8
14.678  36326.7
10      28716.9
6.81292 22417.4
4.64159 17326.6
3.16228 13270.2
2.15443 10111.3
1.4678  7655.31
1       5775.24
0.68129 4316.84
0.46416 3230.06
0.31623 2394.33
0.21544 1769.85
0.14678 1301.58
0.1     950.86
0.06813 688.064
0.04642 494.033
0.03162 352.155
0.02154 249.154
0.01468 176.579
\end{filecontents}


\begin{document}

\begin{tikzpicture}
\begin{loglogaxis}[
    enable tick line clipping=false,
    axis line style=thick,
    width=9cm,
    height=7cm,
    separate axis lines,
    x tick style={black,thick},
    x label style=
        {at={(ticklabel cs:0.5)},anchor=near ticklabel},
    xlabel={$\omega$ [rad/s]},
    xmin=0.01,xmax=100,
    xtick={0.01,0.1,1,10,100},
    xtick pos=bottom,
    xtick align=outside,
    y tick style={black,thick},
    y label style=
        {at={(ticklabel cs:0.5)},anchor=near ticklabel},
    ylabel={$G'$ [Pa]},
    ymin=1,ymax=1000000,
    ytick={1,10,100,1000,10000,100000,1000000},
    ytick pos=left,
    ytick align=outside,
]   

\addplot [only marks,mark=*,mark options={scale=1.3,thick,fill=blue!50}] file {double.dat};

\addplot+[raw gnuplot,red,thick,no marks,smooth] gnuplot {
set log x;
f(x)=\sum_{i=0}^{k} G_k * \frac{x*a_k}{1+(x*a_k)^2};
G=1000;
k=10;
a=1;
fit f(x) 'double.dat' using 1:2 via G_k,a_k,k;
plot [x=0.01:100] f(x);
};


\end{loglogaxis}
\end{tikzpicture}

\end{document}

我可以手动编写 Prony 系列，k=4但我需要使用更多元素，比如说k=100。

Answer 1

对我来说，和的定义不是很清楚，这是一种可能的解释。如果解释是错误的，很容易改变。

\documentclass[tikz,border=3mm]{standalone}
\usepackage{pgfplots}
\usepackage{pgfplotstable}
\usepackage{filecontents}

\pgfkeys{/pgf/number format/.cd,1000 sep={}}




\pgfplotsset{compat=1.17}
\pgfplotsset{label style={font=\Large},
            tick label style={font=\Large}}

\begin{filecontents*}{double.dat}
100     96413.1
68.1292 81638.4
46.4159 68107.2
31.6228 56025.6
21.5443 45408.8
14.678  36326.7
10      28716.9
6.81292 22417.4
4.64159 17326.6
3.16228 13270.2
2.15443 10111.3
1.4678  7655.31
1       5775.24
0.68129 4316.84
0.46416 3230.06
0.31623 2394.33
0.21544 1769.85
0.14678 1301.58
0.1     950.86
0.06813 688.064
0.04642 494.033
0.03162 352.155
0.02154 249.154
0.01468 176.579
\end{filecontents*}

\begin{document}

\begin{tikzpicture}
\pgfplotstableread{double.dat}\mydata
\pgfplotstablegetrowsof{\mydata}%
\pgfmathtruncatemacro{\numrows}{\pgfplotsretval}%
\pgfplotstablegetcolsof{\mydata}%
\pgfmathtruncatemacro{\numcols}{\pgfplotsretval}%
\edef\irun{0}%
\loop
\pgfplotstablegetelem{\irun}{[index]0}\of\mydata%#1=row, #2=column
\edef\currai{\pgfplotsretval}%
\pgfplotstablegetelem{\irun}{[index]1}\of\mydata%#1=row, #2=column
\edef\currGi{\pgfplotsretval}%
\ifnum\irun=0\relax
\edef\PronySum{\currGi*(\currai*x)/(1+\currai*\currai)}%
\else
\edef\PronySum{\PronySum+\currGi*(\currai*x)/(1+pow(\currai*x,2))}%
\fi
\edef\irun{\the\numexpr\irun+1}%
\expandafter\edef\csname PronySum\irun\endcsname{\PronySum}%
\ifnum\irun<\numrows\relax
\repeat
\def\PartialPronySum#1{\csname PronySum#1\endcsname}
\begin{loglogaxis}[%declare function={PronySum(\k)=\PartialPronySum\k;};
    enable tick line clipping=false,
    axis line style=thick,
    width=9cm,
    height=7cm,
    separate axis lines,
    x tick style={black,thick},
    x label style=
        {at={(ticklabel cs:0.5)},anchor=near ticklabel},
    xlabel={$\omega$ [rad/s]},
    xmin=0.01,xmax=100,
    xtick={0.01,0.1,1,10,100},
    xtick pos=bottom,
    xtick align=outside,
    y tick style={black,thick},
    y label style=
        {at={(ticklabel cs:0.5)},anchor=near ticklabel},
    ylabel={$G'$ [Pa]},
    ymin=1,ymax=1000000,
    ytick={1,10,100,1000,10000,100000,1000000},
    ytick pos=left,
    ytick align=outside,
]   

\addplot [only marks,mark=*,mark options={scale=1.3,thick,fill=blue!50}] file {double.dat};
\edef\temp{\noexpand\addplot[domain=0.01:100,red,thick] {\PartialPronySum{10}};}
\temp
\end{loglogaxis}
\end{tikzpicture}

\end{document}