How to integrate one variable in a multivariate formulation+ nonlinear curvefitting
Afficher commentaires plus anciens
Hello everyone,
I am facing a probelm with integration, I want to integrate a multiplication of two terms below is the codes:
function y = MH(para,x)
A = para(1);
mu = para(2);
Ms = para(3);
sigma = para(4);
C = para(5);
kb=1.38065e-23;
T=295;
lognormal =@(d) exp(-(log(d/mu).^2/(2*sigma^2)))./(d*sigma*sqrt((2*pi)));
langevin =@(d,x) coth(((d.^3).*pi/6*Ms*x)/(kb*T))-(kb*T)./(pi/6*Ms*x.*d.^3);
integrand = @(d) ((pi/6)*d.^3).*langevin(d,x).*lognormal(d);
y = A*integral(integrand,0,inf)+C*x;
end
Actually, I am trying to fit a M(H) curve, the formulation is:
where
, 
The unknown parameters are A, Ms, Sigma, Mu and C,.
The two variables are H and d, H is the xdata with a dimension of (61×1) and in this code it is replaced as x, d is the variable that should be integrated and has a limit of integration of (0,inf).
The problem is, the part 'intergrand' isn't integrated at all and it isn't fitted, only the part C*x is fitted, the returned parameter Ms, Sigma and Mu doesn't change from the estimated parameter at beginning, only A and C changed, the result is as below:
Red points are the data that need to be fitted and blue line is the fitting curve.
May I ask what's the reason for that? Thanks for any help in advance!
Réponse acceptée
First, use element-wise operations for everything in the ‘MH’ function.
Second, in the integration use the 'ArrayValued',1 name-value pair.
In ‘integrand’, the argument to lognormal is squared, so I added that, and it seems to work with it, although not without it.
With my synthetic data, I get a decent fit, however the function always returns at least one NaN value. I added the fillmissing call to correct for that in order to test the code, however it may not be necessary with your actual data. If it is, it will be necessary to see what is causing the NaN value. It is not obvious to me what is causing it.
x = (-6:0.24:6)*1E4; % Create Data
y = tanh(x*1E-4) + x*1E-8; % Create Data
B0 = rand(5,1)*100; % Use Appropriate Initial Estimates
B = lsqcurvefit(@MH, B0, x, y)
figure
plot(x, y, '.')
hold on
plot(x, MH(B,x), '-r')
hold off
grid
function y = MH(para,x)
A = para(1);
mu = para(2);
Ms = para(3);
sigma = para(4);
C = para(5);
kb=1.38065e-23;
T=295;
lognormal =@(d) exp(-(log(d/mu).^2./(2*sigma^2)))./(d*sigma*sqrt((2*pi)));
langevin =@(d,x) coth(((d.^3).*pi./6*Ms*x)./(kb*T))-(kb*T)./(pi/6*Ms*x.*d.^3);
integrand = @(d) ((pi/6)*d.^3).*langevin(d,x).*lognormal(d.^2);
y = A*integral(integrand,0,inf, 'ArrayValued',1)+C*x;
y = fillmissing(y,'nearest'); % May Not Be Necessary
end
.
11 commentaires
Hello Sir,
thank you so much for your detailed explanation and effort on solving my problem! I had added the name-value pair on my codes, but it gives a result the same as the former. I tried to fitted the 'lognormal' and 'langevin' seperately, when it comes to 'lognomal', the result is like:
and the value of parameter A, Sigma changes, but I think the integration of the lognormal formulation shouldn't be a line. By 'langevin' only, it can't run and comes erro code says Not enough input arguments. I think may be the expression of the complex formulation has something wrong, but have checked many times.
The full code the program is as below:
xdata = X20210416_Perimag_0_5mg_pro_ml_export(:,1);
ydata = (X20210416_Perimag_0_5mg_pro_ml_export(:,2)-X20210521_Water_1_export(:,2))/1000;
A0 = 1e17;
mu0 = 130e-9;
Ms0 = 3e5;
sigma0 = 0.25;
C0 = 2e-9;
para0 = [A0;mu0;Ms0;sigma0;C0];
para = lsqcurvefit(@func,para0,xdata,ydata);
y_fit = func(para,xdata);
figure;
plot(xdata,ydata,'or');
hold on;
plot(xdata,y_fit,'-b');
function m = func(para,xdata)
m = zeros(size(xdata));
for i = 1:length(xdata)
m(i) = MH(para,xdata(i));
end
end
function y = MH(para,x)
A = para(1);
mu = para(2);
Ms = para(3);
sigma = para(4);
C = para(5);
kb=1.38065e-23;
T=295;
lognormal =@(d) exp(-(log(d/mu).^2./(2*sigma^2)))./(d*sigma*sqrt((2*pi)));
langevin =@(d,x) coth(((d.^3).*pi./6*Ms*x)./(kb*T))-(kb*T)./(pi/6*Ms*x.*d.^3);
integrand = @(d) ((pi/6)*d.^3).*langevin(d,x).*lognormal(d.^2);
y = A*integral(integrand,0,inf,'ArrayValued',1)+C*x;
end
The data for the fitting are enclosures.
If it's possible, could you give me a hand again? Thank you so much for your effort and help!
Hello,
I am not sure about the estimated initial parameter, according to the literature the A should be a value of 10^17, Ms is 10^5, Sigma and C should be the same order as the estimated, Mu is around 10^-9, namely nanoscale. I am sorry that, the value of the data need to be changed as:
xdata = X20210416_Perimag_0_5mg_pro_ml_export(:,1);
ydata = (X20210416_Perimag_0_5mg_pro_ml_export(:,2)-X20210521_Water_1_export(:,2))/1000;
The tested data are without units, so after a convertion I think this order should be right. And after running I get a result as below:
Now the parameter A changes but Ms doesn't, maybe there's still erro by the formulation for fitting?
May I also have a ask, why the quality after fitting isn't so good? Are there other function to get a better result? Thank you so much for your help!
I cannot get it to fit the data at all.
I am not familiar with any of the physics of what you are doing, so I cannot comment on that. Something is obviously wrong somewhere, however I cannot determine what the problem is. The initial parameter estimates could be incorrect, so review those to be certain they are close to what they should be, at least in terms of magnitude.
In any event, I did mange to get your ‘MH’ function to work correctly, and that was what you requested help with.
format short g
Perimag = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/657475/20210416_Perimag_0.5mg%20pro%20ml%20export.txt');
Water = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/657480/20210521%20Water_1%20export.txt');
xdata = Perimag(:,1);
ydata = (Perimag(:,2) - Water(:,2))/1E+3;
% xdata = X20210416_Perimag_0_5mg_pro_ml_export(:,1);
% ydata = (X20210416_Perimag_0_5mg_pro_ml_export(:,2)-X20210521_Water_1_export(:,2))/1000;
A0 = 1e17;
mu0 = 130e-9;
Ms0 = 3e5;
sigma0 = 0.25;
C0 = 2e-9;
para0 = [A0;mu0;Ms0;sigma0;C0];
para = lsqcurvefit(@MH,para0,xdata,ydata)
y_fit = MH(para,xdata);
figure;
plot(xdata,ydata,'or');
hold on;
plot(xdata,y_fit,'-b');
% function m = func(para,xdata)
% m = zeros(size(xdata));
% for i = 1:length(xdata)
% m(i) = MH(para,xdata(i));
% end
% end
function y = MH(para,x)
A = para(1);
mu = para(2);
Ms = para(3);
sigma = para(4);
C = para(5);
kb=1.38065e-23;
T=295;
lognormal =@(d) exp(-(log(d/mu).^2./(2*sigma^2)))./(d*sigma*sqrt((2*pi)));
langevin =@(d,x) coth(((d.^3).*pi./6*Ms*x)./(kb*T))-(kb*T)./(pi/6*Ms*x.*d.^3);
integrand = @(d) ((pi/6)*d.^3).*langevin(d,x).*lognormal(d.^2);
y = A*integral(integrand,0,inf,'ArrayValued',1)+C*x;
end
.
Zhao
le 18 Juin 2021
Hello Sir,
thank you so much for your effort and help! Ya, I need to have a discussion with my superviser about the problem. Actually what I measured is the response of magnetic nanoparticle under magnetic field, so the particle is nanoscale. I had a caculation of the value by hand, A should have a magnitude of 1e17, Mu is median of the partcle and schould be around 130e-9, Sigma should be 1e-2 or 1e-3, Ms should be around 1e5 and C is 1e-9.
Now the problem is, however I change the estimated initial parameters, A, Ms and Sigma don't change at all, even when I get a fitted curve like below:
May I have a ask, what can be the reason of the codes for this problem?
I am sorry that, the litertaure that I reference has failure, it lost a mv=4*pi*10^(-7) in the formulation, the 'langevin' should be:
mv=4*pi*10^(-7);
langevin =@(d,x) coth(((d.^3).*pi/6*Ms*mv*x)/(kb*T))-(kb*T)./(pi/6*Ms*mv*x.*d.^3);
Thank you so much for your help!
Star Strider
le 18 Juin 2021
As always, my pleasure!
Thank you for correcting the ‘langevin’ function and the constant. I will experiment with them later.
The fitted curve appears to be about as good as it can get, with the changes in the ‘langevin’ function and the constant. It will almost never be exact.
If ‘A’, ‘Ms’ and ‘sigma’ do not change in the process of estimating the parameters, it is likely that they are already close to the correct values, or that changes in them do not significantly affect the fit. (I suspect that they do change, although only slightly, so the changes may not be readily apparent.)
.
Star Strider
le 19 Juin 2021
As always, my pleasure!
Thank you for explaining the problems and the approach you are using. If you want tme to experiment further with the fitting, I will. (I do not completely understand what you are doing, since the physics of it are not an area of my expertise.)
Zhao
le 19 Juin 2021
As always, my pleasure!
There must also be other changes that I am not aware of, since even with the ‘langevin’ change, I am not getting a good fit. (I have no idea what ‘mv’ does, or where it should be in the code, so I put it just before ‘langevin’ since that function uses it.)
In any event, parameter estimates and curves fitted to data are approximate, and depend on the data, and to some extent the parameter estimation algorithm itself. Some data willl fit a particular model better than others. (This one appears not to work at all, at least when I try it.)
Please provide the updated and correct function, and other necessary code changes.
format short g
Perimag = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/657475/20210416_Perimag_0.5mg%20pro%20ml%20export.txt');
Water = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/657480/20210521%20Water_1%20export.txt');
xdata = Perimag(:,1);
ydata = (Perimag(:,2) - Water(:,2));
% xdata = X20210416_Perimag_0_5mg_pro_ml_export(:,1);
% ydata = (X20210416_Perimag_0_5mg_pro_ml_export(:,2)-X20210521_Water_1_export(:,2))/1000;
A0 = 1e17;
mu0 = 130e-9;
Ms0 = 3e5;
sigma0 = 0.25;
C0 = 2e-9;
para0 = [A0;mu0;Ms0;sigma0;C0];
para = lsqcurvefit(@MH,para0,xdata,ydata)
y_fit = MH(para,xdata);
figure;
plot(xdata,ydata,'or');
hold on;
plot(xdata,real(y_fit),'-b');
% function m = func(para,xdata)
% m = zeros(size(xdata));
% for i = 1:length(xdata)
% m(i) = MH(para,xdata(i));
% end
% end
function y = MH(para,x)
A = para(1);
mu = para(2);
Ms = para(3);
sigma = para(4);
C = para(5);
kb=1.38065e-23;
T=295;
lognormal =@(d) exp(-(log(d/mu).^2./(2*sigma^2)))./(d*sigma*sqrt((2*pi)));
mv=4*pi*10^(-7);
langevin =@(d,x) coth(((d.^3).*pi/6*Ms*mv*x)/(kb*T))-(kb*T)./(pi/6*Ms*mv*x.*d.^3);
integrand = @(d) ((pi/6)*d.^3).*langevin(d,x).*lognormal(d.^2);
y = A*integral(integrand,0,inf,'ArrayValued',1)+C*x;
end
.
Zhao
le 19 Juin 2021
Star Strider
le 19 Juin 2021
As always, my pleasure!
Please post the code that works!
Plus de réponses (0)
Catégories
En savoir plus sur Noncentral t Distribution dans Centre d'aide et File Exchange
le 17 Juin 2021
le 19 Juin 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
