Approximation of a B-H curve

11 vues (au cours des 30 derniers jours)
Huy
Huy le 14 Déc 2016
Commenté : David Goodmanson le 22 Déc 2016
Hi guys: I'm dealing with an approximation problem. The data of the B-H curve of a core material (steel M4) is given: https://www.dropbox.com/s/r8da63kxg6n9ou4/BH.mat?dl=0
Which then can be plotted using Matlab (blue curve):
I've tried to approximate it with this function (red): B(H)_{app} = 1.62*(1-exp(-H/11.4))
But the result is not good enough for the calculation. The area which is really important to me is the one marked green. I need the approximation in this area to be as good as possible. The slope of the curve is also very important, since dB/dH := the permeability of the material. I'd really appreciate if some one can help me. Thanks a lot.
  2 commentaires
John BG
John BG le 15 Déc 2016
the data in the Dropbox link of your question does not correspond to the data, the original BH graph you show:
load BH.mat
plot(B)
plot(B([2:end]))
B1=1.62*(1-exp(-nB/11.4));hold all;plot(nB,B1);
would you please be so kind to submit the correct data?
awaiting answer
Huy
Huy le 15 Déc 2016
Hi John, thank you so much for responding. I've checked the data and it's actually correct. B values and H vales come in pair and B(H) is a function of H, so you'd have to plot it using plot(H,B) instead of just plot(B); But I've found my answer now thank to David Goodmanson. Thank you anyway for trying to help. Much appreciated.

Connectez-vous pour commenter.

Réponse acceptée

David Goodmanson
David Goodmanson le 15 Déc 2016
Hi Huy, First of all I was able to plot your curve from the .mat file provided, although since it is a short file I think it would have been preferable to just attach it as a text file.
There are a zillion ways to fit a curve like this. One of them is by pade approximation. If you remove the point at the origin from both H and B and fit the rest (point 2 being at around H = 8.1), then
Bfit = polyval(a,H)./polyval(b,H);
with
a =
0.000379890931843
0.087779632450721
-0.400451284821487
0.610758213470415
b =
0.000192370833515
0.054846991096759
-0.055190978462873
1.000000000000000
has a relative error (Bfit-B)/B of less than 3e-5 over the entire range, not counting the point at the origin.
The large number of sig figs in a and b don't really mean anything. The long format was just a convenient way to show their values.
  3 commentaires
Huy
Huy le 21 Déc 2016
Hi David.
I'm so sorry to bother you again, but can you please tell me how you managed to calculate the values of a and b? I've tried polyfit with degree n=3 but somehow get others values than yours. Best regards. Huy
David Goodmanson
David Goodmanson le 22 Déc 2016
Hi Huy, Ordinarily I would say that you could go to Matlab file exchange and pick a pade fit function, but it's not clear if any of them work in the same way as the function below. I have not posted this on the file exchange because I have my doubts whether it is mathematically defensible. On the other hand, it does work a lot of the time. You just have to experiment with na and nb, the degrees of the numerator and denominator polynomials, to see what gives a good result.
function [a,b] = pade(x,y,na,nb)
% Approximate y(x) as a ratio of polynomials in x.
% Num and denom polynomials are of degrees na and nb respectively.
% The approximating function is z = polyval(a,x)./polyval(b,x);
% Coefficient vectors a and b have standard Matlab form
% (highest power first, constant term last), with
% y = Sum{j=1,na+1} aj*x^(na+1-j) / Sum{j=1,nb+1} bj*x^(nb+1-j).
% The constant term in the denominator [i.e. b(end)] equals 1
% by construction. If nb is not specified, then nb = na.
% a and b have the same form (row or column) as y.
%
% [a,b] = pade(x,y,na,nb)
% David Goodmanson
% for email address: a = char(abs('NR[\Rij1kT\de%[hg')+(22:-1:6))
isrowy = isrow(y);
x = x(:);
y = y(:);
if nargin == 3, nb = na; end
% create vandermonde matrices va and -y.* vb
v = [ones(size(x)) cumprod(repmat(x,1,max(na,nb)),2)];
v = [v(:,1:na+1), repmat(-y,1,nb).*v(:,2:nb+1)];
% solve for a,b coefficients
ab = v\y;
a = flipud(ab(1:na+1));
b = flipud([1; ab(na+2:end)]);
if isrowy
a = a.'; % row vectors
b = b.';
end

Connectez-vous pour commenter.

Plus de réponses (1)

Image Analyst
Image Analyst le 15 Déc 2016
Then why don't you just pass it the data in the green zone? Don't pass it ALL the data, just pass it ONLY the data you care about.

Catégories

En savoir plus sur Get Started with Curve Fitting Toolbox dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by