I want to find the value of two constants from an equation with two variables

I am trying to find the values of the constant elastic modulus (E1) and viscosity modulus (eta1) from the Kelvin creep recovery equation. The data points of strain and time needs to be exported from an excel file(given) and the value of constant stress is given as 6 nN.

The equation being used is: here sigma(0) is the constant stress(given) and epsilon is the variable strain(data points given)

Answer 1

[Star Strider]

@NDC —

To estimate the initial parameters for this, since is known, the asymptote as will be , and that gives a way to estimate E. Here, is essentially a time constant, so if you can calculate the time at 0.6 times the asymptote (use interp1), you can estimate the initial value of η as well, having already calcualted E.

Use those as initial values with the code you already wrote.

.

[Nidhi]

Thank you. I will try doing that.

[Star Strider]

As always, my pleasure!

.

[Nidhi]

@Star Strider Hi. I tried solving the equation using this approach. Is it correct or am I solving it wrongly?

clc
clear
v = xlsread('Book1.xlsx');
x=v(:,2);
y=v(:,1);
%to find initial parameters taking into account two initial data points of
%the variables
syms E n
eqns=[(6/E)*(1-exp((-E.*v(1,1))/n))==v(1,2),(6/E)*(1-exp((-E.*v(2,1))/n))==v(2,2)];
S=vpasolve(eqns,[E n]);
S.E
S.n
%to find value of constants taking into account all the data points for the
%variables
x = @(p,v) (6/p(1))*(1-exp(-(p(1).*v(:,1))/p(2)));
p0= [0.03934;0.000068];    %initial parameters                                                                
SSECF = @(p) sum((v(:,2) - x(p,v)).^2);                                             
[abc, SSE] = fminsearch(SSECF, p0);                                              
E = abc(1)
n = abc(2)

[Star Strider]

‘Is it correct or am I solving it wrongly?’

Not quite, and not quite.

I made a few corrections to the initial parameter estimate calculations, and slighty changed the objective function (added an offset, since the dependent variable does not start at 0) —

v = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/707317/Book1.xlsx');
x=v(:,2);
y=v(:,1);

% figure
% plot(y,x,'.')
% grid
% xlabel('y')
% ylabel('x')
% text(y,x,compose('%02d',1:numel(x)))
%to find initial parameters taking into account two initial data points of
%the variables
syms E n
eqns=[(6/E)*(1-exp((-E.*v(1,1))/n))==v(1,2),(6/E)*(1-exp((-E.*v(2,1))/n))==v(2,2)];
S=vpasolve(eqns,[E n]);
S.E

ans = 

0.039344262295081967213114751056344

S.n

ans = 

0.000068043909801043736024547852887911

%to find value of constants taking into account all the data points for the
%variables
xfcn = @(p,v) (6/p(1)).*(1-exp(-(p(1).*v)./p(2))) + p(3);
% p0= [0.03934;0.000068];    %initial parameters          

idxrng = find(x < max(x)*0.999);
tau = interp1(x(idxrng)-min(x), y(idxrng), 0.6*(max(x)-min(x)));
sigma0 = 6;
p0(1) = sigma0/x(end);                                          % Estimate 'E' (Asymptote = 'x(end)')
p0(2) = tau;                                                    % Extimate 'eta'
p0(3) = min(x)                                                  % Offset

p0 = 1×3

0.0382 0.7724 152.0300

SSECF = @(p) sum((x - xfcn(p,y)).^2);                                             
[abc, SSE] = fminsearch(SSECF, p0)                                              

abc = 1×3

1.1248 0.8704 151.8299

SSE = 0.1148

E = abc(1)

E = 1.1248

n = abc(2)

n = 0.8704

figure
plot(y, x, 'p')
hold on
plot(y, xfcn(abc,y), '-r')
hold off
grid
text(2.5, 154.5, sprintf('$\\epsilon(t) = \\frac{%.0f}{%.2f}(1-e^{\\frac{%.2f}{%.2f} t})$',6,E,E,n), 'Interpreter','latex', 'FontSize',14)

That appears to me to be a decent fit.

.

[Nidhi]

@Star Strider Thank you so much! Can you please explain the use of this two lines? Also is 0.6 given because the value of sigma0 is 6?

idxrng = find(x < max(x)*0.999);
tau = interp1(x(idxrng)-min(x), y(idxrng), 0.6*(max(x)-min(x)));

[Star Strider]

Sure!

The ‘idxrng’ variable creates a vector of subscripts that select a part of the signal vector that do not contain any duplicates in the first argument. This is necessary in order for interp1 to work. It interpolates to find the approximate time constant, τ.

The actual value of the time constant, τ, is :

tau = 1-exp(-1)

tau = 0.6321

I refer you to the Wikipedia article on Time constant for a full explanation.

The value of 0.6 is close enough for our purposes, since we are using it as an initial parameter estimate.

.

[Nidhi]

@Star Strider Thank you!

[Star Strider]

As always, my pleasure!

.

Answer 2

[darova]

Try fit

[Nidhi]

@darova Hi. I have never previously used "fit" and am quite new to MATLAB. Can you kindly, please show me how I can implement "fit" in this case.