0

I am trying to fit Prony series to my viscoelasticity data in PGFPlots. The problem I am running with is how to mathematically define the Prony series that gnuplot can recognise the function and do the fitting. The Prony series is \sum_{i=0}^{k} G_k * \frac{x*a_k}{1+(x*a_k)^2}. Anyone know how to do this? Here is the complete code I am using:

\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}

I can manually write the Prony series for k=4 but I need to use more elements, say k=100.

  • Please make the question a bit clearer. What are the G_k and a_k? The k-th entries in your table? And i seems to be your summation index. However, the argument of the sum does not depend on i. – user194703 Apr 16 at 6:00
  • The function definition in the gnuplot program is no valid gnuplot code, i expect sum [<var> = <start> : <end>] <expression>! Do you want to use latex to generate the Prony series for use in the gnuplot program? Please minimize your example to the core problem as far as possible. – hefe Apr 16 at 6:30
  • @Schrödinger's cat Terribly sorry for the mistake. It must be G_i and a_i, these are the parameters for generalised Maxwell model. Also, i=0 cannot be, must start from 1. The correct function to be fitted is \sum_{i=1}^{k} \frac{G_i*x*a_i}{1+(x*a_i)^2}. – Gagik Apr 16 at 19:54
  • Well, the answer below was assuming this. (Note that pgfplots starts numbering its rows at 0.) – user194703 Apr 16 at 20:29
  • @hefe I used the following expression f(x)=sum[i=1:100] G_i*x*a_i/(1+(x*a_i)^2) you suggested. Texpad compiles it without errors. The only trouble is that the gnuplot doesn't take the summation over x=1:100 range. It only calculates for one element x=1. Any thoughts on that? – Gagik Apr 19 at 0:01
0

To me the definition of the sum is not very clear, this is a possible interpretation. If the interpretation is wrong, it is easy to change this.

\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}

enter image description here

| improve this answer | |
  • I don't see the fitting part in this answer. However, no meaningful response is possible until the ambiguities in the question is resolved. – hefe Apr 16 at 8:04
  • @hefe This answer enables the OP to compute sums over expressions which draw the coefficients from a table. No more, no less. If they come up with a formula that provides one with a proper fit, then this answer can be used to visualize the result. To me it is not even clear if this it supposed to be a fitting function. – user194703 Apr 16 at 8:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.