Need help to fit the data without error.
2 次查看(过去 30 天)
显示 更早的评论
% 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).');
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.
4 个评论
回答(1 个)
Torsten
2024-3-4
编辑:Torsten
2024-3-5
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)
norm(fitfun(initialGuess,rData,tData))
norm(fitfun(fit,rData,tData))
[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
9 个评论
Torsten
2024-3-13
编辑:Torsten
2024-3-13
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Get Started with Curve Fitting Toolbox 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!