Approximation of a B-H curve

Question

Huy 2016-12-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/316928-approximation-of-a-b-h-curve

评论： David Goodmanson 2016-12-22

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 个评论
显示无隐藏无

John BG 2016-12-15

在 MATLAB Online 中打开

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

John BG <mailto:jgb2012@sky.com jgb2012@sky.com>

Huy 2016-12-15

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.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

David Goodmanson 2016-12-15

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/316928-approximation-of-a-b-h-curve#answer_247258

在 MATLAB Online 中打开

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 个评论
显示 1更早的评论隐藏 1更早的评论

Huy 2016-12-21

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 2016-12-22

在 MATLAB Online 中打开

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