How to curve fit an equation that gets values from another equation?

7 views (last 30 days)
I have 38 data points for x and t (input variables) and 38 data points for y (output variable). In equation 1, a and b are coefficients to be found. z is the output variable which iterates and gets input from its previous value. That means when i=1, z_2 gets value from z_1. when i=2, z_3 gets value from z_2, and so on.
We can take z_0 = 0.5
.................equation (1)
solving equation 1 should give me 38 points for z, which will depend on the x and t.
....................equation (2)
In equation 2, x is the input data, y(fit) is the output and z comes from equation 1. c,d and e are coefficient to be found.
Now, I want to fit equation 2 to y vs x curve (which I have from experiments) and find the values of the constants a,b,c,d and e.
I have previously used lsqcurvefit to fit a single equation. But in this case there are 2 equations and the second one gets values from the first one.
Any help with this will be much appreciated. Thanks in advance.
  5 Comments
Ases Akas Mishra
Ases Akas Mishra on 17 Feb 2022
Hi, I am new to MATLAB central and I dont know how to load the data. I have pasted the data below. Let me know if you need it in some other format.
x = 1, 1.44, 2.07, 2.98, 4.28, 6.16, 8.86, 12.7, 18.3, 26.4, 37.9, 54.6, 78.5, 113, 162, 234, 336, 483, 695, 1000, 1000, 695, 483, 336, 234, 162, 113, 78.5, 54.6, 37.9, 26.4, 18.3, 12.7, 8.86, 6.16, 4.28, 2.98, 2.07
t = 1800, 3600, 5400, 7200, 9000, 10800, 12600, 14400, 16200, 18000, 19800, 21600, 23400, 25200, 27000, 28800, 30600, 32400, 34200, 36000, 37800, 39600, 41400, 43200, 45000, 46800, 48600, 50400, 52200, 54000, 55800, 57600, 59400, 61200, 63000, 64800, 66600, 68400
y = 1.419, 1.3369, 1.1017, 1.4513, 2.2703, 3.1121, 4.1855, 4.2835, 4.6603, 4.8933, 5.287, 5.877, 6.4844, 7.2159, 7.9723, 12.472, 15.202, 21.45, 24.348, 27.307, 26.128, 21.671, 18.912, 16.709, 14.279, 12.437, 10.866,9.5601, 8.5242, 7.733, 7.0755, 6.5612, 6.1244, 5.531, 4.8106, 3.2184, 2.1077, 1.6007
And here is my code-
clc
clear all
close all
x= xlsread('data.xlsx','0.5','B4:B41');
y= xlsread('data.xlsx','0.5','C4:C41');
t = xlsread('data.xlsx','0.5','H4:H41');
for i = 1 : length(x)
% to find z
a = 100; % I have used some random value for a and b. But they should be obtained from the fit.
b = 0.05;
A(i,:) = a/((b.*x(i)) + a); % parts of equation 1
e(i,:) = -t(i) / (a.*x(i)+b);
z(1,1) = 0.5; %initialize z
z(i+1,:) = A(i) + (z(i) - A(i)) * exp(e(i)); %equation 1
z1{i,:} = z(i+1,:) %to remove the extra term(z=0.5)
z_final = cell2mat(z1); %converting to double
% equation 2
c = 0.5; %Again guessing values for c,d,e
d = 2.2;
e = 0.35;
y_fit(i,:) = (c + d.*z_final(i)).*x(i)^e %equation 2
end

Sign in to comment.

Answers (3)

Rik
Rik on 18 Feb 2022
There is probably something wrong with your loop, as you're using e both in the first part as an array, and in the second part. To avoid errors I renamed the one in the first part to E. Since you're only ever accessing e(i) anyway, there is also no need for indexing (same for A).
This is all assuming your function actually works as expected. Tweak your function as needed.
x = [1, 1.44, 2.07, 2.98, 4.28, 6.16, 8.86, 12.7, 18.3, 26.4, 37.9, 54.6, 78.5, 113, 162, 234, 336, 483, 695, 1000, 1000, 695, 483, 336, 234, 162, 113, 78.5, 54.6, 37.9, 26.4, 18.3, 12.7, 8.86, 6.16, 4.28, 2.98, 2.07];
t = [1800, 3600, 5400, 7200, 9000, 10800, 12600, 14400, 16200, 18000, 19800, 21600, 23400, 25200, 27000, 28800, 30600, 32400, 34200, 36000, 37800, 39600, 41400, 43200, 45000, 46800, 48600, 50400, 52200, 54000, 55800, 57600, 59400, 61200, 63000, 64800, 66600, 68400];
y = [1.419, 1.3369, 1.1017, 1.4513, 2.2703, 3.1121, 4.1855, 4.2835, 4.6603, 4.8933, 5.287, 5.877, 6.4844, 7.2159, 7.9723, 12.472, 15.202, 21.45, 24.348, 27.307, 26.128, 21.671, 18.912, 16.709, 14.279, 12.437, 10.866,9.5601, 8.5242, 7.733, 7.0755, 6.5612, 6.1244, 5.531, 4.8106, 3.2184, 2.1077, 1.6007];
%initial estimates
a = 100;
b = 0.05;
c = 0.5;
d = 2.2;
e = 0.35;
lsqcurvefit(@(abcde,x)find_y_fit(x,t,abcde(1),abcde(2),abcde(3),abcde(4),abcde(5)),...
[a b c d e],x,y)
Solver stopped prematurely. lsqcurvefit stopped because it exceeded the function evaluation limit, options.MaxFunctionEvaluations = 5.000000e+02.
ans = 1×5
83.0941 -0.0005 -337.2435 339.7435 0.2726
function y_fit=find_y_fit(x,t,a,b,c,d,e)
%it is possible to optimize/clean this a lot futher, but let's first get it working
y_fit=zeros(size(x));
z=zeros(size(x));
z(1,1) = 0.5; %initialize z
for n = 1 : numel(x)
% to find z
A = a/((b.*x(n)) + a); % parts of equation 1
E = -t(n) / (a.*x(n)+b);
z(n+1,:) = A + (z(n) - A) * exp(E); %equation 1
z1{n,:} = z(n+1,:); %to remove the extra term(z=0.5)
z_final = cell2mat(z1); %converting to double
% equation 2
y_fit(n) = (c + d.*z_final(n)).*x(n)^e ; %equation 2
end
end
  3 Comments
Rik
Rik on 20 Feb 2022
I have little experience with advanced curve fitting. My naive approach would be to use something like ndgrid to generate a few initial estimates for the 5 parameters. Note that this will blow up quickly, as you will have to run the fit for N^5 iterations for N points per parameter.
Since you didn't provide the code you used, I can't tell whether the points simply lack sorting, or if there is something else wrong. You have not confirmed my edits to your code are actually correct. You really should. Use the debugger to step through the code line by line.

Sign in to comment.


Mathieu NOE
Mathieu NOE on 18 Feb 2022
hello
my first job was to make your code a bit more efficient and streamlined
also notice you use both "e" for variable name and for a constant , avoid this as it can lead to errors and unpredictable results. I renamed the "e" from your equation 1 as ""ee" to avoid any conflict with the constant "e".
so bellow is basically your code a bit rewanped :
clc
clear all
close all
% x= xlsread('data.xlsx','0.5','B4:B41');
% y= xlsread('data.xlsx','0.5','C4:C41');
% t = xlsread('data.xlsx','0.5','H4:H41');
data= xlsread('Data.xlsx');
x = data(:,1);
t = data(:,2);
y = data(:,3);
z = 0.5;
a = 100; % I have used some random value for a and b. But they should be obtained from the fit.
b = 0.05;
c = 0.5; %Again guessing values for c,d,e
d = 2.2;
e = 0.35;
for i = 1 : length(x)
A(i,:) = a/((b.*x(i)) + a); % parts of equation 1
ee(i,:) = -t(i) / (a.*x(i)+b);
z(i+1,:) = A(i) + (z(i) - A(i)) * exp(ee(i)); %equation 1
end
z_final = z(2:end); %to remove the (first) extra term(z=0.5)
% equation 2
y_fit = (c + d.*z_final).*(x.^e); %equation 2
plot(t,y,t,y_fit);
second step, your equations are now a function to minimize the distance (norm) vs the experimental data , and you get the optimal parameters in array opt_vector ( that contains your five constants [a b c d e] ):
here I used the fminsearch function, but you can use more advanced optmization function if you need to add some constraints on the parameters (I don't have the Optimization Toolbox);
results :
opt_vector = 12.3055 -0.0016 -0.3599 1.9294 0.3889
code :
clc
clear all
close all
% x= xlsread('data.xlsx','0.5','B4:B41');
% y= xlsread('data.xlsx','0.5','C4:C41');
% t = xlsread('data.xlsx','0.5','H4:H41');
data= xlsread('Data.xlsx');
x = data(:,1);
t = data(:,2);
y = data(:,3);
%% fit equation to data
% parameters
a = 100; % I have used some random value for a and b. But they should be obtained from the fit.
b = 0.05;
c = 0.5; %Again guessing values for c,d,e
d = 2.2;
e = 0.35;
params = [a b c d e];
opt_vector = fminsearch(@ (params) norm (y - myequation(params, data)) ,params);
y_fitted = myequation(opt_vector, data);
figure(1);
plot(t,y,t,y_fitted);
legend('raw','fit');
function y_fit = myequation(params, data)
x = data(:,1);
t = data(:,2);
y = data(:,3);
% parameters
a = params(1);
b = params(2);
c = params(3);
d = params(4);
e = params(5);
% init z
z = 0.5;
for i = 1 : length(x)
A(i,:) = a/((b.*x(i)) + a); % parts of equation 1
ee(i,:) = -t(i) / (a.*x(i)+b);
z(i+1,:) = A(i) + (z(i) - A(i)) * exp(ee(i)); %equation 1
end
z_final = z(2:end); %to remove the (first) extra term(z=0.5)
% equation 2
y_fit = (c + d.*z_final).*(x.^e); %equation 2
end
  3 Comments
Ases Akas Mishra
Ases Akas Mishra on 20 Feb 2022
Follow up request- I would also like to collect the final values of z_final and y_fit (after optimization) so that I can use it for other purposes such as plotting using another software.

Sign in to comment.


Alex Sha
Alex Sha on 18 Feb 2022
How about the results below:
Residual Sum of Squares (RSS): 34.9122402573179
Root of Mean Square Error (RMSE): 0.958510910040288
Correlation Coef. (R): 0.991557637426582
R-Square: 0.983186548338985
Parameter Best Estimate
-------------------- -------------
a 4.63375351494044E-8
b -1.56073777954123E-7
c 1.80838380873546
d -1.35268042423128
e 0.387487908425502
  1 Comment
Ases Akas Mishra
Ases Akas Mishra on 20 Feb 2022
Hi ALex
This looks fine. Like I have mentioned before, the issue is to somehow guess the initial values of a,b,c,d,e such that the fit is perfect. There will be an unique combination of a,b,c,d,e when the fit will be perfect and then the 'calculated y' curve won't be symmetric (along vertical axis), like it is now. Please share your code. Thanks for your interest.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!