Using least square method , determining unknown values ?

Hello, I have a question which i couldn't solve. Please help me.

The table shows the experimental results of the measured Coulomb Force, F, between two charges, (q1 and q2) corresponding to distance r. General form of the Coulomb Force is:

 F = (1/4*pi*ε0)*q1*q2/(r^n)

Write a program to determine values of ε0 and n using least-square method.

 r(cm) | F (N) 
-------+------------ 
 40    |  2.5 +- 0.1 
 30    |  4.5 +- 0.1 
 20    | 10.1 +- 0.2 
 10    | 40.5 +- 0.4 
  5    |162.0 +- 0.9

q1 = 5 μC q2 = 9 μC

[Wayne King]

What have you done to solve your homework problem?

Show the MATLAB code you have written and any error messages you are getting.

[osman]

I used curve fitting tool to determine n and epsilon and i did but i couldn't write a program determining it. So i can't show any program . Any clue about writing this program would be appreciated.

Answer 1

[Friedrich]

Hi,

is this even doable with normal least square? This is a nonlinear problem. Normally you would try to write the problem as

Ax = b

and than do

(A'A)^-1 A'b

But in your case you can't write the problem in that form since you have it like this

(1/4*pi*q1*q2) * (ε0/r^n) = F

Where r changed from equation to equation. You would need to use the logarithm to get the n but than the ε0 gets into a logarithm statement.

You would need to get it in a form which has a form like this

A * [ε0;n] = F

(sorry for the bad formatting, but writing formulas is pretty bad on answers). Where A is a matrix and F a vector.

I think one can't split ε0 and n up in such a way. So this seem that some algorithm for a non linear problem must be used here.

[osman]

ε0 and n actually are not changing. n is 2 and ε0 is constant value. But i should show they are.

[Friedrich]

I never stated that n and ε0 change. But with a small but not very good program (used optimization toolbox) I get ε0 = 114.6834 n = 2.0005 as coeffs. The result strongly depends on the start values.

I created a function which depends on ε0 and n in which I calculate ||f(ε0,n) -F||.

function err = my_func( x )

%x(1) = ε0

%x(2) = n;

fix = 1/4*pi*5*9*x(1); %1/4*pi*q1*q2*ε0

%calc f - F

f_F = [ fix/40^x(2) - 2.5

fix/30^x(2) - 4.5

fix/20^x(2) - 10.1

fix/10^x(2) - 40.5

fix/5^x(2) - 162];

err = norm(f_F);

end

Than I used fminunc to calculate the minimum:

%guess 150 for ε0 and 1 for n

[coeefs,fval] = fminunc(@my_func,[150,1])

%test if the values are okay

disp([num2str(1/4*pi*5*9*coeefs(1)/40^coeefs(2)) ' = 2.5'])

disp([num2str(1/4*pi*5*9*coeefs(1)/30^coeefs(2)) ' = 4.5'])

disp([num2str(1/4*pi*5*9*coeefs(1)/20^coeefs(2)) ' = 10.1'])

disp([num2str(1/4*pi*5*9*coeefs(1)/10^coeefs(2)) ' = 40.5'])

disp([num2str(1/4*pi*5*9*coeefs(1)/5^coeefs(2)) ' = 162'])

Which leads to a warning of the solver and these results:

coeefs =

114.6834 2.0005

fval =

0.0368

2.5287 = 2.5

4.4961 = 4.5

10.1182 = 10.1

40.4867 = 40.5

162.0019 = 162

[Friedrich]

if you have the global optimization toolbox you can use the follow code where the startpoint doesn't matter:

rs = RandomStartPointSet('NumStartPoints',1000); %generate 1000 startpoints

options = optimset('TolFun',1e-16,'TolX',1e-16,'MaxFunEvals',1000,'MaxIter',1000);

problem = createOptimProblem('fminunc','objective',@my_func,'x0',[1,100],'options',options) %choose x0 which is far from the actual solution

[xmin,fmin,flag,outpt,allmins] = run(MultiStart,problem,rs)