5
$\begingroup$

I'm trying to fit experimental data with FindFit and implicitly calculated function. It seems, that function works fine and satisfactorily approximates data when parameters are fed manually. However when trying to apply FindFit with the same set of parameters as starting values, it returns "Tensors {..1..} and {..2..} have incompatible shapes". These tensors are equal. Manual calculation shows, that their elements are differences between experimental and fitted values (for starting values of parameters) for each experimental point. Thus by multiplying them, FindFit probably tries to calculate the sum of squared residuals. If above is correct, why does this error emerge? I tried substituting the fitting function with simple second-order polynomial and it works. Here is the code:

Remove["Global`*"]    
Hk = 952.45844; Rahe = 815; const = 797; AngCr = 
     54.35 \[Degree]; Ang = 45 \[Degree];
    a[b_?NumericQ] := ArcSin[Tan[b]/Tan[AngCr] ];
    NewH[b_?NumericQ, c_?NumericQ, Ang_?NumericQ, Hamp_?NumericQ] := 
      Flatten[({
          {Cos[c], -Sin[c], 0},
          {Sin[c], Cos[c], 0},
          {0, 0, 1}
         }).({
          {Cos[b], 0, Sin[b]},
          {0, 1, 0},
          {-Sin[b], 0, Cos[b]}
         }).({
          {1, 0, 0},
          {0, Cos[a[b]], -Sin[a[b]]},
          {0, Sin[a[b]], Cos[a[b]]}
         }).({
          {Hamp*Cos[Ang \[Degree]]},
          {Hamp*(-Sin[Ang \[Degree]])},
          {0}
         })];
    (* Implicit fitting function: *)
    Rfit[b_?NumericQ, c_?NumericQ, Ang_?NumericQ, Hamp_?NumericQ, 
      Rphe_?NumericQ] := 
     Rahe*Sin[theta] - 
       Rphe*Cos[theta]^2*
        Abs[(NewH[b, c, Ang, Hamp][[1]]*NewH[b, c, Ang, Hamp][[2]])/(
         NewH[b, c, Ang, Hamp][[1]]^2 + NewH[b, c, Ang, Hamp][[2]]^2)] + 
       const /. 
      Quiet@NSolve[
        Tan[theta] == (Hk*Sin[theta] + NewH[b, c, Ang, Hamp][[3]])/Sqrt[
          NewH[b, c, Ang, Hamp][[1]]^2 + NewH[b, c, Ang, Hamp][[2]]^2] && 
         theta <= Pi/2 && theta > 0.5, theta, Reals]
    FindFit[ExpData, 
     Rfit[b, c, Ang, Hamp, 
      Rphe], {{b, 0, -0.1, 0.1}, {c, 0, -0.1, 0.1}, {Rphe, 450, 400, 
       500}}, Hamp]

    (* Manual approximation works fine: *)
    Show[
     ListPlot[
      Parallelize[
       Table[Flatten[{Hamp, 
          Rfit[-2.5 \[Degree], 0 \[Degree], 45, Hamp, 465]}], {Hamp, -401,
          401, 10}]], PlotStyle -> {Red, AbsolutePointSize[6]}],
     ListPlot[ExpData, PlotStyle -> {Black, AbsolutePointSize[5]}],
     PlotRange -> All
     ]

Here is the piece of experimental data:

ExpData = {{-377.4285`, 1503.685`}, {-366.8757`, 
    1509.651`}, {-357.7828`, 1514.686`}, {-347.0045`, 
    1520.59`}, {-338.1022`, 1525.638`}, {-327.2154`, 
    1530.05`}, {-318.4717`, 1534.859`}, {-307.7064`, 
    1539.4`}, {-299.1506`, 1543.492`}, {-288.0877`, 
    1548.19`}, {-279.8359`, 1552.482`}, {-268.8763`, 
    1557.043`}, {-260.4716`, 1560.633`}, {-249.6292`, 
    1563.926`}, {-240.9715`, 1567.46`}, {-230.3227`, 
    1570.869`}, {-220.9978`, 1573.732`}, {-209.8838`, 
    1578.208`}, {-202.0967`, 1580.598`}, {-191.8868`, 
    1583.704`}, {-181.9674`, 1586.1`}, {-172.2505`, 
    1588.278`}, {-162.1441`, 1591.064`}, {-152.7795`, 
    1592.871`}, {-142.4023`, 1594.542`}, {-133.2651`, 
    1597.427`}, {-125.6026`, 1597.731`}, {-113.7561`, 
    1601.733`}, {-103.0388`, 1603.187`}, {-94.00935`, 
    1604.345`}, {-83.34793`, 1605.78`}, {-74.33023`, 
    1606.449`}, {-63.44335`, 1608.238`}, {-54.81877`, 
    1608.54`}, {-44.04633`, 1608.988`}, {-35.41862`, 
    1610.753`}, {-24.58675`, 1610.793`}, {-17.20935`, 
    1610.65`}, {-5.281529`, 1612.059`}, {2.261039`, 
    1610.745`}, {13.69646`, 1611.337`}, {23.07399`, 
    1611.651`}, {33.15822`, 1609.854`}, {43.12982`, 
    1609.793`}, {52.70468`, 1608.815`}, {62.87678`, 
    1608.523`}, {72.25062`, 1607.767`}, {82.43324`, 
    1606.551`}, {91.69639`, 1605.46`}, {97.61601`, 
    1601.687`}, {111.4225`, 1602.904`}, {121.9307`, 
    1601.799`}, {130.9983`, 1599.582`}, {141.4585`, 
    1598.096`}, {150.7068`, 1595.881`}, {161.252`, 
    1593.875`}, {170.2011`, 1591.611`}, {180.8177`, 
    1588.666`}, {189.7459`, 1585.693`}, {200.3986`, 
    1583.918`}, {209.3968`, 1581.302`}, {219.8005`, 
    1578.068`}, {227.8856`, 1574.321`}, {239.7572`, 
    1571.973`}, {249.5857`, 1568.514`}, {259.3328`, 
    1565.229`}, {269.4772`, 1561.71`}, {279.0915`, 
    1558.152`}, {289.2914`, 1553.891`}, {298.7476`, 
    1549.908`}, {308.9017`, 1545.625`}, {318.3074`, 
    1541.339`}, {328.4131`, 1536.248`}, {337.7596`, 
    1532.195`}, {347.9528`, 1528.05`}, {357.3661`, 
    1522.416`}, {367.6034`, 1517.481`}, {377.0643`, 1512.4`}};

Thanks a lot in advance.

Improve this question
$\endgroup$
4
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – Dr. belisarius
    Commented Apr 1, 2015 at 12:52
  • $\begingroup$ Note, that you require solutions that are Reals. If for your parameter set or during the fitting procedure no real solution is found, this will probably result in a problem down the road. This might cause your incompatible tensors. $\endgroup$
    – mikuszefski
    Commented Apr 1, 2015 at 15:51
  • $\begingroup$ Also note that you give Ang in degrees already, but later define a function using Ang Degree. This is probably not wanted and in general dangerous causing many errors difficult to find. Try not to use a degree conversion inside your functions, but only on function call, i.e. on parameter definition. $\endgroup$
    – mikuszefski
    Commented Apr 1, 2015 at 15:56
  • $\begingroup$ @mikuszefski thanks for the Degree comment, missed it, now it's fixed. As for Reals, I believe it's OK, because I tuned the Rfit function so that it returns exactly one real value (for sets of parameters, which are experimentally reasonable). $\endgroup$
    – Alexander
    Commented Apr 2, 2015 at 9:45

1 Answer 1

Reset to default
11
$\begingroup$

Some errors corrected and a few tricks

  • Using memoization,
  • "Precomputing" the solutions by using Reduce once instead of NSolve each time,
  • Rewriting the functions in a compact and more understandable way,
  • Reducing the overall calcs needed,
  • Using NonlinearModelFit with Method -> {"NMinimize", Method -> "NelderMead"} for efficiency
  • Reduced the ConfidenceLevel for an additional speedup

: 

Hk = 952.45844; Rahe = 815; const = 797; AngCr = 54.35 °; Ang = 45 °;
a[b_] := ArcSin[Tan[b]/Tan[AngCr]];
red = Quiet@Reduce[Tan[t] == (Hk*Sin[t] + bb)/cc && cc > 0, t, Reals] /. C[1] -> 0;
rt[b_, c_] := First@Select[t/. {red/. {bb-> b, cc-> c}//N//ToRules}, 1/2<#<Pi/2 &]

rid = Reverse@IdentityMatrix@3;
NewH[b_, c_, w_, Hamp_] := NewH[b, c, w, Hamp] =
  Hamp Evaluate[Dot@@MapThread[RotationMatrix, {{c, b, a@b}, rid}].{Cos@w,-Sin@w,0}]

Rfit[b_?NumericQ, c_, Ang_, Hamp_?NumericQ, Rphe_] := Rfit[b, c, Ang, Hamp, Rphe] =
  Module[{v = NewH[b, c, Ang, Hamp], sq, t},
   sq = v[[1]]^2 + v[[2]]^2;
   t = rt[v[[3]], Sqrt@sq];
   Rahe*Sin[t] - Rphe*Cos[t]^2*Abs[v[[1]] v[[2]]/sq] + const]

Now:

k = NonlinearModelFit[ExpData, {Rfit[b, c, Ang, Hamp, Rphe], 
                     -.1 < b < .1 && -.1 < c < .1 && 400 < Rphe < 500}, 
                       {b, c, Rphe}, Hamp, 
                       Method -> {"NMinimize", Method -> "NelderMead"}, 
                       ConfidenceLevel -> .5]

A fair fit:

ListLinePlot@k[{"Response", "PredictedResponse"}]

Mathematica graphics

The parameters:

k["BestFitParameters"]
(* {b -> -0.0471051, c -> 0.1, Rphe -> 499.547}*)

You may (or may not) get better results by increasing the confidence level at the expense of more calc time.

Improve this answer
$\endgroup$
17
  • $\begingroup$ To simplify even further, note the NewH[] function is basically a coordinate transformation of a vector. You can remove the Flatten[] and the inner braces in the last vector. $\endgroup$
    – mikuszefski
    Commented Apr 1, 2015 at 16:00
  • $\begingroup$ The main error, I guess was the fact that the function returned {value}, which you corrected by First[]. With this adjustment FindFit[] should work as well, right? Moreover, this makes it a error due to erroneous function definition and not a FindFit problem. This is already suggested by the error message when running the original post. Finally, looking at the boundary conditions of @Alexander, one also has to notice that c and Rphe are on the limit here, rendering the result questionable (Not from the answer point of view, but from the physics point of view...just saying). $\endgroup$
    – mikuszefski
    Commented Apr 1, 2015 at 17:33
  • $\begingroup$ @mikuszefski the main issue was calc time. I massaged the fitting procedure until it spitted out something in a few minutes. $\endgroup$
    – Dr. belisarius
    Commented Apr 1, 2015 at 18:12
  • $\begingroup$ @mikuszefski Regarding the "edge values" for the parameters, I really don't know. Without a deeper understanding of the physics it's impossible to say if those intervals are right $\endgroup$
    – Dr. belisarius
    Commented Apr 1, 2015 at 18:14
  • $\begingroup$ Nice answer & speed-up, +11 $\endgroup$
    – ciao
    Commented Apr 1, 2015 at 23:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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