I would like to reproduce the results from a publication (https://doi.org/10.1007/s00161-015-0415-8) for verification. For this I want to solve the following system of differential equations:
$$\dot{x} = \frac{\beta\left(...\right)}{x_0}\left(\mu_{ma} - \left(1-x_0\right)\mu_{ra} - x_0\mu_{c}-\frac{R\Theta}{M}\ln{\frac{x}{x_0-x}}\right),$$ $$\dot{\delta} = -A\left(...\right)\left(1-\frac{x}{x_0}\right) \left(d_{ma}\delta - w_{ma}p + e_{ma}\left(\Theta-\Theta_0\right)\right)$$ with $$\beta = \biggl|\frac{x-x_{eq}}{x_{eq}}\biggl|\beta_0\exp{\frac{c_1\left(\Theta-\Theta_0\right)}{c_2+\Theta-\Theta_0}},$$ $$A = A_0\exp{\frac{c_1\left(\Theta-\Theta_0\right)}{c_2+\Theta-\Theta_0}}$$ and $$\mu_{ma}\left(p, \Theta, \delta\right) = \mu_{ma0} + d_{ma}\frac{\delta^2}{2} - \frac{p^2}{2\rho_{ma}\kappa_{ma}} - s_{ma0} \left(\Theta-\Theta_0\right) + e_{ma} \left(\Theta-\Theta_0\right)\delta + \frac{\alpha_{ma}}{\rho_{ma}}p\left(\Theta-\Theta_0\right) - w_{ma}p\delta - \frac{\beta_{ma0}}{2}\left(\Theta-\Theta_0\right)^2 - \left(c_{ma}-\beta_{ma0}\Theta_0\right) \left(\Theta\ln{\frac{\Theta}{\Theta_0}}-\left(\Theta-\Theta_0\right)\right),$$
$$\mu_{ra}\left(p, \Theta, \delta\right) = \mu_{ra0} - s_{ra0} \left(\Theta-\Theta_0\right) - \frac{\beta_{ra0}}{2}\left(\Theta-\Theta_0\right)^2 - \left(c_{ra}-\beta_{ra0}\Theta_0\right) \left(\Theta\ln{\frac{\Theta}{\Theta_0}}-\left(\Theta-\Theta_0\right)\right),$$
$$\mu_{c}\left(p, \Theta, \delta\right) = \mu_{c0} - \frac{p^2}{2\rho_{c}\kappa_{c}} - s_{c0} \left(\Theta-\Theta_0\right) + \frac{\alpha_{c}}{\rho_{c}}p\left(\Theta-\Theta_0\right) - \left(c_{c}-\beta_{c}\Theta_0\right) \left(\Theta\ln{\frac{\Theta}{\Theta_0}}-\left(\Theta-\Theta_0\right)\right) - \frac{\beta_{c}}{2}\left(\Theta-\Theta_0\right)^2.$$
Values for the parameters can be taken from the following picture:
The following equations with the initial values $\Theta=500 K$ and $p=0$ are used as initial values:
$$\delta_{eq} = \frac{1}{d_{ma}}\left(w_{ma}p-e_{ma}\left(\Theta-\Theta_0\right)\right),$$ $$x_{eq} = \frac{x_0}{1+\exp{\left(-\frac{M\left(\mu_{ma}\left(p,\Theta,\delta_{eq}\right)-x_0\mu_c\left(p,\Theta\right)-\left(1-x_0\right)\mu_{ra}\left(\Theta\right)\right)}{R\Theta}\right)}}.$$
As well as the following curves for $\Theta$ as a function of time t:
Curves for theta as a function of time
I use in Matlab the solver ODE15s with the following settings:
temp_rate1 = 0.01;
temp_rate2 = 0.01;
Tmax = 500;
Tmin=300 ;
y0=[x_init delta_init];
zyk=2;
ta=0;
te=0.5*zyk*(Tmax-Tmin)*(1/abs(temp_rate1)+1/abs(temp_rate2));
dt=(te/zyk)/80000;
tspan = ta:dt:te;
options=odeset('RelTol',1e-12,'AbsTol',[1e-12 1e-12],'InitialStep',dt/3,'MaxStep',dt);
[t,y] = ode15s('DGL',tspan,y0, options);
Unfortunately, I do not get the results shown in the publication. Below you can see the results from the publication and my generated results. The solver does not converges at a certain point. Furthermore, the typical plateau formation is not reached in the solution I generated. In contrast to the results from the publication, my results show qualitatively the same course for each cooling rate.
own result of x for Coolingrate 0.01 K/s
own result of x for Coolingrate 0.1 K/s
own result of x for Coolingrate 1 K/s
Can anyone help me further in solving this DGL system?
Below you will find excerpts from my Matlab code for the DGL.
Thanks a lot!
Felix
t_mid = 200/Cooling_rate
if t < t_mid
Theta = Tmax - Cooling_rate*t;
else
Theta = Tmin - Cooling_rate*t_mid + Heating_rate*t;
end
delta_eq = 1/d_ma*(w_ma*p-e_ma*(Theta-Theta_0));
mu_ma = mu_ma0 + d_ma*(y(2)^2)/2 - ...
(p^2)/(2*rho_ma*kappa_ma) - s_ma0*(Theta-Theta_0)+...
e_ma*(Theta-Theta_0)*y(2) + (alpha_ma/rho_ma)*p*...
(Theta-Theta_0)-w_ma*p*y(2) - (beta_ma/2)*...
(Theta-Theta_0)^2-(c_ma-beta_ma*Theta_0)*...
(Theta*log(Theta/Theta_0)-(Theta-Theta_0));
mu_ra = mu_ra0- s_ra0*(Theta-Theta_0)-(beta_ra/2)*...
(Theta-Theta_0)^2 -(c_ra-beta_ra*Theta_0)*...
(Theta*log(Theta/Theta_0)-(Theta-Theta_0));
mu_c = mu_c0 - (p^2)/(2*rho_c*kappa_c) - s_c0*(Theta-Theta_0) + ...
(alpha_c/rho_c)*p*(Theta-Theta_0) - (c_c-beta_c*Theta_0)*...
(Theta*log(Theta/Theta_0)-(Theta-Theta_0))-beta_c/2*(Theta-Theta_0)^2;
mu_ma_delta_eq = mu_ma0 + d_ma*(delta_eq^2)/2 - ...
(p^2)/(2*rho_ma*kappa_ma) - s_ma0*(Theta-Theta_0)+...
e_ma*(Theta-Theta_0)*delta_eq + (alpha_ma/rho_ma)*p*...
(Theta-Theta_0)-w_ma*p*delta_eq - (beta_ma/2)*...
(Theta-Theta_0)^2-(c_ma-beta_ma*Theta_0)*...
(Theta*log(Theta/Theta_0)-(Theta-Theta_0));
x_eq = x_0/(1+exp(-(M*(mu_ma_delta_eq-x_0*mu_c-...
(1-x_0)*mu_ra))/(R*Theta)));
A = A_0*exp((c_1*(Theta-Theta_0))/(c_2+Theta-Theta_0));
beta = abs((Theta-x_eq)/x_eq)^Chi*...
beta_0*exp((c_1*(Theta-Theta_0))/(c_2+Theta-Theta_0));
% x_dot
x_dot = (beta/x_0)*(mu_ma-(1-x_0)*mu_ra-x_0*mu_c-...
((R*Theta)/M)*log(y(1)/(x_0-y(1))));
%delta_dot
delta_dot = -A*(1-(y(1)/x_0))*(d_ma*y(2)-w_ma*p+e_ma*(Theta-Theta_0));
%function for ODE-Solver
dydt = zeros(2,1);
dydt(1) = x_dot;
dydt(2) = delta_dot;
you may check your implementation of $\beta$ in the matlab code. Original source states $x-x_\text{eq}$. However you evaluating $\theta - x_\text{eq}$. Maybe you will be able to get the desired solution. If not, you may contact the authors. As I understand, the corresponding author is professor at an university in Munich. Usually, they will provide you with their source code.