1
$\begingroup$

I have written the following code in Mathematica

f1[u_?NumberQ, r_?NumberQ] := 
  Y /. 
    FindRoot[
      1 - (1 - u) Y - u (1 + 3.14^2/(r Y^((2 0.63 - 1)/0.62))^2)^(1/2) Y^(0.63/.52), 
      {Y, 0.000001, 1}];

rf1[u_?NumberQ, r_?NumberQ] := 
  3.14^2/(Sqrt[r f1[u, r]^((2 0.63 - 1)/0.62)]^4) 
    (1 + 3.14^2/Sqrt[r f1[u, r]^((2 0.63 - 1)/0.62)]^2)^(-1) 
    (0.63/0.52 + ((1 - u) f1[u, r])/(1 - (1 - u) f1[u, r]))^(-1);

rf2[u_?NumberQ, r_?NumberQ] := -((2 0.63 - 1)/0.52) r f1[u, r]^((2 0.63 - 1)/0.62);

f2[u_?NumberQ, r_?NumberQ] := 
  1/0.52*(rf1[u, r] r f1[u, r]^((2 0.63 - 1)/0.62))/(1 + rf1[u, r] rf2[u, r])

w0 = -1;
w1 = 1.5 ;
w2 = -1.4;
w3 = -0.8;
w4 = 0.07;
b = 2.691047;

w[theta_?NumberQ, l0_?NumberQ, m0_?NumberQ] := 
  m0 l0 (w0 + w1 theta^2 + w2 theta^4 + w3 theta^6 + w4 theta^8);
wp[theta_, l0_, m0_] := D[w[theta, l0, m0], theta];

l[theta_?NumericQ, l0_?NumericQ] := l0 theta (1 - theta^2);
lp[theta_?NumericQ, l0_?NumericQ] := -2 l0 theta^2 + l0 (1 - theta^2);
k[theta_?NumericQ] := (1 - b theta^2);
kp[theta_?NumericQ] := -3.382094 theta;

data = 
  {{127481.8183296798`, 0.001`}, {1990.5248645895465`, 0.011`}, 
   {650.2050275373107`, 0.021`}, {331.742608489351`, 0.031`}, 
   {204.75942189403582`, 0.041`}, {140.54519314852138`, 0.051000000000000004`}, 
   {103.24315592880238`, 0.061`}, {79.5018435395663`, 0.07100000000000001`}, 
   {63.37868794896939`, 0.081`}, {51.88495362079122`, 0.091`}}

m[theta_?NumericQ, r_?NumericQ, u_?NumericQ, l0_?NumericQ, m0_?NumericQ] := 
  ((2 - 0.11 f2[u, r]) w[theta, l0, m0] kp[theta] - 
   wp[theta, l0, m0] k[theta])/(lp[theta, l0] k[theta] - 
   (3/2 + 0.325 4.8 - 3/2) f2[u, r] l[theta, l0] kp[theta]);
cdd1[theta_, r_, u_, l0_, m0_] := Block[{x}, ND[m[theta, x, u, l0, m0], x, r]]

At this point when I use

FindFit[data, cdd1[1, r, u , l0, m0] , {u, l0, m0}, r]

But I get various errors for example "1 is not a valid variable.", it seems as I am not experienced in Mathematica, even with NonlinearModelFit I get error, I don't know where I made mistake. Could you please help me? Thanks

$\endgroup$
1
$\begingroup$

The primary error here is from your definition of wp (and also through w). You can't take the derivative with respect to $\theta$ having put $\theta=1$ into the equation already. Try replacing them with:

w[theta_, l0_?NumberQ, m0_?NumberQ] := 
  m0 l0 (w0 + w1 theta^2 + w2 theta^4 + w3 theta^6 + w4 theta^8);

wp[theta_, l0_, m0_] := D[w[th, l0, m0], th] /. th -> theta;

(There are other ways to do this for wp).

This lets me evaluate cdd1 for numerical values successfully.

$\endgroup$
0
$\begingroup$

these defs should not be be SetDelayed or NumericQ :

w[theta_, l0_, m0_] = 
  m0 l0 (w0 + w1 theta^2 + w2 theta^4 + w3 theta^6 + w4 theta^8);
wp[theta_, l0_, m0_] = D[w[theta, l0, m0], theta];

then you can not (I think) use FindFit where the fit parameters need to be numeric. Try this:

FindMinimum[
 Total[(cdd1[1, #[[1]], u, l0, m0] - #[[2]])^2 & /@ data ] , {{u, 
   1}, {l0, 1}, {m0, 1}}]

{0.00512945, {u -> 0.819404, l0 -> 13062.9, m0 -> 164.599}}

Edit: yes you can use FindFit. It takes a lot longer to give a similarly poor result.

FindFit[data, cdd1[1, x, u, l0, m0], {{u, 1}, {l0, 1}, {m0, 1}}, x]

the fit quality here is not good, but it may be your function just doesn't fit, or maybe you need to find better start parameters.

Edit-the whole code

Needs["NumericalCalculus`"] (*for ND*)
f1[u_?NumberQ, r_?NumberQ] := 
  Y /. FindRoot[
    1 - (1 - u) Y - 
     u (1 + 3.14^2/(r Y^((2 0.63 - 1)/0.62))^2)^(1/
         2) Y^(0.63/.52), {Y, 0.000001, 1}];
rf1[u_?NumberQ, r_?NumberQ] := 
  3.14^2/(Sqrt[r f1[u, r]^((2 0.63 - 1)/0.62)]^4) (1 + 
      3.14^2/Sqrt[r f1[u, r]^((2 0.63 - 1)/0.62)]^2)^(-1) (0.63/
       0.52 + ((1 - u) f1[u, r])/(1 - (1 - u) f1[u, r]))^(-1);
rf2[u_?NumberQ, 
   r_?NumberQ] := -((2 0.63 - 1)/0.52) r f1[u, r]^((2 0.63 - 1)/0.62);
f2[u_?NumberQ, r_?NumberQ] := 
 1/0.52*(rf1[u, r] r f1[u, r]^((2 0.63 - 1)/0.62))/(1 + 
     rf1[u, r] rf2[u, r])
w0 = -1;
w1 = 1.5;
w2 = -1.4;
w3 = -0.8;
w4 = 0.07;
b = 2.691047;
w[theta_, l0_, m0_] = 
  m0 l0 (w0 + w1 theta^2 + w2 theta^4 + w3 theta^6 + w4 theta^8);
wp[theta_, l0_, m0_] = D[w[theta, l0, m0], theta];
l[theta_, l0_] = l0 theta (1 - theta^2);
lp[theta_, l0_] = -2 l0 theta^2 + l0 (1 - theta^2);
k[theta_] = (1 - b theta^2);
kp[theta_] = -3.382094 theta;
data = {{127481.8183296798`, 0.001`}, {1990.5248645895465`, 
   0.011`}, {650.2050275373107`, 0.021`}, {331.742608489351`, 
   0.031`}, {204.75942189403582`, 0.041`}, {140.54519314852138`, 
   0.051000000000000004`}, {103.24315592880238`, 
   0.061`}, {79.5018435395663`, 
   0.07100000000000001`}, {63.37868794896939`, 
   0.081`}, {51.88495362079122`, 0.091`}}
m[theta_?NumericQ, r_?NumericQ, u_?NumericQ, l0_?NumericQ, 
   m0_?NumericQ] := ((2 - 0.11 f2[u, r]) w[theta, l0, m0] kp[theta] - 
     wp[theta, l0, m0] k[theta])/(lp[theta, l0] k[
       theta] - (3/2 + 0.325 4.8 - 3/2) f2[u, r] l[theta, l0] kp[
       theta]);
cdd1[theta_, r_, u_, l0_, m0_] := 
 Block[{x}, ND[m[theta, x, u, l0, m0], x, r]]
fit = FindMinimum[
  Total[(cdd1[1, #[[1]], u, l0, m0] - #[[2]])^2 & /@ data], {{u, 
    1}, {l0, 1}, {m0, 1}}]
$\endgroup$
  • $\begingroup$ @ george2079 Thanks, but I got this problem FindMinimum::nrlnum: The function value {1.41421 (-12.+GlobalND[m[1.,x,1.,1.,1.],x,1.]),1.41421 (-10.+GlobalND[m[1.,x,1.,1.,1.],x,1.9]),1.41421 (-8.2+GlobalND[m[1.,x,1.,1.,1.],x,2.6]),1.41421 (-6.9+GlobalND[m[1.,x,1.,1.,1.],x,3.4]),1.41421 (-5.9+Global`ND[m[1.,x,1.,1.,1.],x,5.])} is not a list of real numbers with dimensions {5} at {u,l0,m0} = {1.,1.,1.}. >> could youplease help me? $\endgroup$ – SF.AF Dec 13 '17 at 16:43
  • $\begingroup$ i pasted the whole code I ran into the answer. maybe you need to restart your kernel. $\endgroup$ – george2079 Dec 13 '17 at 16:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.