Need help to fit the data without error.

% Data
rData = [5.3571429, 0.096854535; 10.714286, 0.055104186; 16.071429, 0.042811499; 21.428571, 0.024825886; 26.785714, 0.023279183; 32.142857, 0.016328542; 37.5, 0.0092185037; 42.857143, 0.0075624777; 48.214286, 0.0023514323; 53.571429, 0.001637045; 58.928571, -0.0024887011; 64.285714, -0.0034741333; 69.642857, -0.0056340032; 75, -0.0040906991; 80.357143, -0.0039738424; 85.714286, -0.0044593789; 91.071429, -0.0054884315; 96.428571, -0.0037277341; 101.78571, -0.0041691748; 107.14286, -0.0039292558; 112.5, -0.0037408923; 117.85714, -0.0040700255; 123.21429, -0.0028904555; 128.57143, -0.0022557232; 133.92857, -0.0020756487; 139.28571, -0.0020739949; 144.64286, -0.0015149035; 150, -0.0019796368; 155.35714, -0.00068430865; 160.71429, -0.00060721168; 166.07143, -0.00055972397; 171.42857, -0.0011788755; 176.78571, -0.00090675531; 182.14286, -0.00060012026; 187.5, 7.6071311e-6];
tData = [5.3571429, 0.081473653; 10.714286, -0.0076210718; 16.071429, -0.038565046; 21.428571, -0.014000405; 26.785714, -0.042161254; 32.142857, -0.071404281; 37.5, -0.066992712; 42.857143, -0.031355057; 48.214286, -0.02043848; 53.571429, -0.025259291; 58.928571, -0.019615094; 64.285714, -0.015185751; 69.642857, -0.012213914; 75, -0.0047624032; 80.357143, -0.00041652762; 85.714286, 0.0028162852; 91.071429, 0.00979253; 96.428571, 0.0080315783; 101.78571, 0.0034739882; 107.14286, 0.0021786814; 112.5, 0.0043349925; 117.85714, 0.0053397331; 123.21429, 0.0061087654; 128.57143, 0.0028425693; 133.92857, 0.002129577; 139.28571, 0.0068534431; 144.64286, 0.0071201038; 150, 0.0099290536; 155.35714, 0.0089545127; 160.71429, 0.0079282308; 166.07143, 0.0075533041; 171.42857, 0.01092774; 176.78571, 0.012219652; 182.14286, 0.01013098; 187.5, 0.0096622622];
% Define equations
syms x;
mu = sym('mu'); lambda = sym('lambda'); ke = sym('ke'); ko = sym('ko');
Alpha = @(mu, lambda) 1/(2*mu + lambda);
r = @(x) x;
theta = @(x) x;
% Define equations
eqn1 = x^2*diff(r(x), x, x) + x*diff(r(x), x) - r(x) + Alpha(mu, lambda)^-1*ke^2*x^2*r(x) - Alpha(mu, lambda)^-1*ko^2*x^2*theta(x) == 0;
eqn2 = x^2*diff(theta(x), x, x) + x*diff(theta(x), x) - theta(x) + mu^-1*ko^2*x^2*r(x) + mu^-1*ke^2*x^2*theta(x) == 0;
% Convert symbolic equations to function handles
eqn1_func = matlabFunction(eqn1);
eqn2_func = matlabFunction(eqn2);
% Define the model
funr = @(params, x) deval(ode45(@(x, y) [eqn1_func(params(1), params(2), params(3), params(4), x); eqn2_func(params(1), params(2), params(3), params(4), x)], [0.75, 187.5], [0.1625, 0]), x, 1);
% Fit the data with initial guess values
initialGuess = [15, 50, 0.01, 0.01];
fit = lsqcurvefit(@(params, x) funr(params, x), initialGuess, rData(:,1).', rData(:,2).');

Error using symengine>@(ke,ko,mu,x)(ke.^2.*x.^3)./mu+(ko.^2.*x.^3)./mu==0.0
Too many input arguments.

Error in solution>@(x,y)[eqn1_func(params(1),params(2),params(3),params(4),x);eqn2_func(params(1),params(2),params(3),params(4),x)] (line 17)
funr = @(params, x) deval(ode45(@(x, y) [eqn1_func(params(1), params(2), params(3), params(4), x); eqn2_func(params(1), params(2), params(3), params(4), x)], [0.75, 187.5], [0.1625, 0]), x, 1);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in solution>@(params,x)deval(ode45(@(x,y)[eqn1_func(params(1),params(2),params(3),params(4),x);eqn2_func(params(1),params(2),params(3),params(4),x)],[0.75,187.5],[0.1625,0]),x,1) (line 17)
funr = @(params, x) deval(ode45(@(x, y) [eqn1_func(params(1), params(2), params(3), params(4), x); eqn2_func(params(1), params(2), params(3), params(4), x)], [0.75, 187.5], [0.1625, 0]), x, 1);

Error in solution>@(params,x)funr(params,x) (line 20)
fit = lsqcurvefit(@(params, x) funr(params, x), initialGuess, rData(:,1).', rData(:,2).');

Error in lsqcurvefit (line 256)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

fit2 = lsqcurvefit(@(params, x) funr(params, x), initialGuess, tData(:,1).', tData(:,2).');
% Plot the fitted functions
x_values = linspace(1, 187.5, 1000);
r_fit = funr(fit, x_values);
theta_fit = funr(fit2, x_values);
figure;
plot(x_values, r_fit, 'r', 'LineWidth', 2);
hold on;
scatter(rData(:,1), rData(:,2), 'b');
xlabel('x');
ylabel('r(x)');
title('Fitted r(x) vs. Data');
legend('Fitted r(x)', 'Data');
figure;
plot(x_values, theta_fit, 'g', 'LineWidth', 2);
hold on;
scatter(tData(:,1), tData(:,2), 'b');
xlabel('x');
ylabel('theta(x)');
title('Fitted theta(x) vs. Data');
legend('Fitted theta(x)', 'Data');

Getting folowwing errors:

Error using symengine>@(ke,ko,mu,x)(ke.^2.*x.^3)./mu+(ko.^2.*x.^3)./mu==0.0

Too many input arguments.

Error in test5>@(x,y)[eqn1_func(params(1),params(2),params(3),params(4),x);eqn2_func(params(1),params(2),params(3),params(4),x)] (line 21)

funr = @(params, x) deval(ode45(@(x, y) [eqn1_func(params(1), params(2), params(3), params(4), x); eqn2_func(params(1), params(2), params(3), params(4), x)], [0.75, 187.5], [0.1625, 0]), x, 1);

Error in odearguments (line 92)

f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in test5>@(params,x)deval(ode45(@(x,y)[eqn1_func(params(1),params(2),params(3),params(4),x);eqn2_func(params(1),params(2),params(3),params(4),x)],[0.75,187.5],[0.1625,0]),x,1) (line 21)

funr = @(params, x) deval(ode45(@(x, y) [eqn1_func(params(1), params(2), params(3), params(4), x); eqn2_func(params(1), params(2), params(3), params(4), x)], [0.75, 187.5], [0.1625, 0]), x, 1);

Error in test5>@(params,x)funr(params,x) (line 25)

fit = lsqcurvefit(@(params, x) funr(params, x), initialGuess, rData(:,1).', rData(:,2).');

Error in lsqcurvefit (line 225)

initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Error in test5 (line 25)

fit = lsqcurvefit(@(params, x) funr(params, x), initialGuess, rData(:,1).', rData(:,2).');

Caused by:

Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

[Torsten]

Look at the examples here

https://uk.mathworks.com/help/symbolic/odetovectorfield.html

to avoid the errors you made above.

[Star Strider]

Note that ‘eqn1’ and ‘eqn2’ are not differential equations —

% Define equations
syms x;
mu = sym('mu'); lambda = sym('lambda'); ke = sym('ke'); ko = sym('ko');
Alpha = @(mu, lambda) 1/(2*mu + lambda);
r = @(x) x;
theta = @(x) x;
% Define equations
eqn1 = x^2*diff(r(x), x, x) + x*diff(r(x), x) - r(x) + Alpha(mu, lambda)^-1*ke^2*x^2*r(x) - Alpha(mu, lambda)^-1*ko^2*x^2*theta(x) == 0

eqn1 = 

eqn2 = x^2*diff(theta(x), x, x) + x*diff(theta(x), x) - theta(x) + mu^-1*ko^2*x^2*r(x) + mu^-1*ke^2*x^2*theta(x) == 0

eqn2 = 

What do you want to do with them, and what parameters do you want to estimate?

.

[tuhin]

But eqn 1 and 2 are differential equations. This is the actual equations:

(*Define the equations*)

\[Alpha][\[Mu]_, \[Lambda]_] := 1/(2*\[Mu] + \[Lambda]);

eqns[\[Mu]_, \[Lambda]_, ke_, ko_] := {

x^2*r''[x] + x*r'[x] - r[x] + \[Alpha][\[Mu], \[Lambda]]^-1*ke^2*x^2*r[x] - \[Alpha][\[Mu], \[Lambda]]^-1*ko^2*x^2*\[Theta][x] == 0,

x^2*\[Theta]''[x] + x*\[Theta]'[x] - \[Theta][x] + \[Mu]^-1*ko^2*x^2*r[x] + \[Mu]^-1*ke^2*x^2*\[Theta][x] == 0};

I want to use and solve this two coupled differential eqns (with four boundary conditions) and fit it with the experimental data of r(x) vs x and \theta(x) vs x. For this I want to tune those four parameters mu; lambda; ke; ko. Please let me know if you need more clarifications.

[Walter Roberson]

Your latex doesn't make sense.

Answer 1

[Torsten]

I assumed you start integration at x = 5.3571429.

You will need to specify dr/dx and dtheta/dx in the y0-vector below.

rData = [5.3571429, 0.096854535; 10.714286, 0.055104186; 16.071429, 0.042811499; 21.428571, 0.024825886; 26.785714, 0.023279183; 32.142857, 0.016328542; 37.5, 0.0092185037; 42.857143, 0.0075624777; 48.214286, 0.0023514323; 53.571429, 0.001637045; 58.928571, -0.0024887011; 64.285714, -0.0034741333; 69.642857, -0.0056340032; 75, -0.0040906991; 80.357143, -0.0039738424; 85.714286, -0.0044593789; 91.071429, -0.0054884315; 96.428571, -0.0037277341; 101.78571, -0.0041691748; 107.14286, -0.0039292558; 112.5, -0.0037408923; 117.85714, -0.0040700255; 123.21429, -0.0028904555; 128.57143, -0.0022557232; 133.92857, -0.0020756487; 139.28571, -0.0020739949; 144.64286, -0.0015149035; 150, -0.0019796368; 155.35714, -0.00068430865; 160.71429, -0.00060721168; 166.07143, -0.00055972397; 171.42857, -0.0011788755; 176.78571, -0.00090675531; 182.14286, -0.00060012026; 187.5, 7.6071311e-6];
tData = [5.3571429, 0.081473653; 10.714286, -0.0076210718; 16.071429, -0.038565046; 21.428571, -0.014000405; 26.785714, -0.042161254; 32.142857, -0.071404281; 37.5, -0.066992712; 42.857143, -0.031355057; 48.214286, -0.02043848; 53.571429, -0.025259291; 58.928571, -0.019615094; 64.285714, -0.015185751; 69.642857, -0.012213914; 75, -0.0047624032; 80.357143, -0.00041652762; 85.714286, 0.0028162852; 91.071429, 0.00979253; 96.428571, 0.0080315783; 101.78571, 0.0034739882; 107.14286, 0.0021786814; 112.5, 0.0043349925; 117.85714, 0.0053397331; 123.21429, 0.0061087654; 128.57143, 0.0028425693; 133.92857, 0.002129577; 139.28571, 0.0068534431; 144.64286, 0.0071201038; 150, 0.0099290536; 155.35714, 0.0089545127; 160.71429, 0.0079282308; 166.07143, 0.0075533041; 171.42857, 0.01092774; 176.78571, 0.012219652; 182.14286, 0.01013098; 187.5, 0.0096622622];
initialGuess = [15, 50, 0.01, 0.01,...
                rData(1,2),(rData(2,2)-rData(1,2))/(rData(2,1)-rData(1,1)),...
                tData(1,2),(tData(2,2)-tData(1,2))/(tData(2,1)-tData(1,1))];
options = optimset('MaxFunEvals',10000,'MaxIter',10000);
fit = lsqnonlin( @(params)fitfun(params,rData,tData), initialGuess,[0 0 0 0 -Inf -Inf -Inf -Inf],[],options)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

fit = 1×8

1.0e+03 * 2.6450 0.0499 0.0000 0.0000 0.0001 -0.0000 0.0000 -0.0000

norm(fitfun(initialGuess,rData,tData))

ans = 3.4186

norm(fitfun(fit,rData,tData))

ans = 0.1558

[T,Y] = odefun(fit,rData,tData);
figure(1)
hold on
plot(rData(:,1),rData(:,2),'o')
plot(T,Y(:,1))
hold off

figure(2)
hold on
plot(tData(:,1),tData(:,2),'o')
plot(T,Y(:,3))
hold off

function res = fitfun(params,rData,tData)
  [T,Y] = odefun(params,rData,tData);
  res = [rData(:,2)-Y(:,1);tData(:,2)-Y(:,3)];
end

function [T,Y] = odefun(params,rData,tData)
  mu = params(1);
  lambda = params(2);
  ke = params(3);
  ko = params(4);
  fun = @(x,y)[y(2);(-x*y(2)+y(1)-ke^2*x^2/(2*mu+lambda)*y(1)+ko^2*x^2*y(3)/(2*mu+lambda))/x^2;...
               y(4);(-x*y(4)+y(3)-ko^2/mu*x^2*y(1)-ke^2/mu*x^2*y(3))/x^2];
  tspan = rData(:,1);
  y0 = [params(5);params(6);params(7);params(8)]; %[r,dr/dx,theta,dtheta/dx]
  [T,Y] = ode45(fun,tspan,y0);
end

[tuhin]

There is no integration, only differentiation. Can you please write full code here so that I just can check?

[Torsten]

There is no integration, only differentiation.

?? You have differential equations for r and theta. Thus you will have to integrate to compare with rData and tData .

Can you please write full code here so that I just can check?

Shouldn't be too much work for you to supply dr/dx and dtheta/dx in the y0 vector :-)

[tuhin]

Thanks, Torsten for your valuable suggestions!

[tuhin]

Getting same parameter values from those two fittings.

Fitting values for the first fitting:

mu = 0.51423

lambda = 0.00019766

ke = 0.030559

ko = 2.3197e-14

-----------------------------------

Fitting values for the second fitting:

mu = 0.51423

lambda = 0.00019766

ke = 0.030559

ko = 2.3197e-14. I do not think that is right. Can you please look into this?

[Torsten]

I changed

fun = @(x,y)[y(2);(-x*y(2)+y(1)-(2*mu+lambda)*ke^2*x^2*y(1)+(2*mu+lambda)*ko^2*x^2*y(3))/x^2;...
             y(4);(-x*y(4)+y(3)-ko^2/mu*x^2*y(1)-ke^2/mu*x^2*y(3))/x^2];

to

fun = @(x,y)[y(2);(-x*y(2)+y(1)-ke^2*x^2/(2*mu+lambda)*y(1)+ko^2*x^2*y(3)/(2*mu+lambda))/x^2;...
             y(4);(-x*y(4)+y(3)-ko^2/mu*x^2*y(1)-ke^2/mu*x^2*y(3))/x^2];

Because you didn't supply an initial value vector y0 for the equations, I defined them also as parameters to be fitted. If you know initial values, you should change the y0 vector accordingly.

If you are not satisfied with the fit or the values of the parameters, you should play with the initialGuess vector or/and set bounds on the parameters in the arrays lb and ub in the call to lsqnonlin.

Good luck.

[tuhin]

I have a doubt. I am getting two different sets of Ke and Ko values for r and theta fitting plots. If I want to get the same set of Ke and Ko values for the two fittings (r and theta), how should I modify it? Any suggestion from you Torsten?

[Torsten]

If you have, e.g. two, different sets of r and t data and you get different fit parameters for these data sets if you fit them separately, you should make one big dataset of the two and fit them simultaneously. Or do you mean something else ?

[tuhin]

Till now, I was varying those four parameters of mu, lamda, ke and ko for these two fittings (i.e. r(x) vs x and t(x) vs x). That means I am getting Ke_r, Ko_r, mu_r and lamda_r for r(x) vs x. For, t(x) vs x I am getting ke_t, ko_t, mu_t and lamda_t. Now I want the ke_r should be same as ke_t and ko_r should be same as ko_t. For that I want to change the script.

[Torsten]

You have one fit in your code above from which you get Ke, Ko, mu and lamda by solving 4 differential equations for r, r', t and t' simultaneously. I don't know what you mean by

Till now, I was varying those four parameters of mu, lamda, ke and ko for these two fittings (i.e. r(x) vs x and t(x) vs x). That means I am getting Ke_r, Ko_r, mu_r and lamda_r for r(x) vs x. For, t(x) vs x I am getting ke_t, ko_t, mu_t and lamda_t.