Message Boards Message Boards

0
|
5237 Views
|
6 Replies
|
1 Total Likes
View groups...
Share
Share this post:

[?] Minimize a very complicated expression?

Hello! I am trying to make a new model which will fit some experimental data. My analytic formulas depend on two parameters :a and b, and I want to find the values for a and b such that the expressions will fit the experimental data in the best way possible.

You can see that the expressions I used are quite complicated and it would take too much space to write them as a single term, so I decided to write multiple terms and the final one is expressed in terms of the others (which is a correct approach I suppose).

For finding the parameters I opted for the chi^2 function minimum, which you can compute easily with "Minimize" . Unfortunately, when I try to run it, it doesn't work: I get only the expression of chi^2 in terms of a and b without any numerical values.

What could be the problem?

Than you in advance! I have also added the .nb file for accessibility.

Data sets

In[1]:= tsd1 = {1739.3, 1936, 2199.3, 2514.1, 2900.4, 3350.7, 3866, 4444.6, 5083.6, 
   5780.6, 6533.2, 7338.8, 8196.5, 9106.2, 10068.7, 11085.2, 12156.3, 13282.5,
    14461.8, 15688.8, 16957.8, 18261.3};
tsd1 = tsd1/1000;
tsd2 = {3078.9, 3486.3, 3957.9, 4492.2, 5088, 5742.6, 6453.8, 7220, 8039.9, 
   8912.8, 9839.3, 10819.5, 11854.2, 12943.1, 14086.1, 15283.4, 16530.9};
tsd2 = tsd2/1000;
tsd3 = {3863.3, 4368.8, 4936.8, 5563.6, 6248.7, 6989.9, 7785.8, 8635.6, 
   9538.1, 10493.9, 11503.1, 12566.1, 13678.5, 14825.8};
tsd3 = tsd3/1000;
tsd4 = {6319.5, 6964.5, 7666.7, 8421.3, 9231.4, 10096.7, 11017.2, 11992.9, 
   13024.5, 14110};
tsd4 = tsd4/1000;

In[9]:= data1 = {{6.5, 1.7393}, {8.5, 1.936}, {10.5, 2.1993}, {12.5, 2.5141}, {14.5, 
    2.9004}, {16.5, 3.3507}, {18.5, 3.866}, {20.5, 4.4446}, {22.5, 
    5.0836}, {24.5, 5.7806}, {26.5, 6.5332}, {28.5, 7.3388}, {30.5, 
    8.1965}, {32.5, 9.1062}, {34.5, 10.0687}, {36.5, 11.0852}, {38.5, 
    12.1563}, {40.5, 13.2825}, {42.5, 14.4618}, {44.5, 15.6888}, {46.5, 
    16.9578}, {48.5, 18.2613}};
data2 = {{13.5, 3.0789}, {15.5, 3.4863}, {17.5, 3.9579}, {19.5, 
    4.4922}, {21.5, 5.088}, {23.5, 5.7426}, {25.5, 6.4538}, {27.5, 
    7.22}, {29.5, 8.0399}, {31.5, 8.9128}, {33.5, 9.8393}, {35.5, 
    10.8195}, {37.5, 11.8542}, {39.5, 12.9431}, {41.5, 14.0861}, {43.5, 
    15.2834}, {45.5, 16.5309}};
data3 = {{16.5, 3.8633}, {18.5, 4.3688}, {20.5, 4.9368}, {22.5, 
    5.5636}, {24.5, 6.2487}, {26.5, 6.9899}, {28.5, 7.7858}, {30.5, 
    8.6356}, {32.5, 9.5381}, {34.5, 10.4939}, {36.5, 11.5031}, {38.5, 
    12.5661}, {40.5, 13.6785}, {42.5, 14.8258}};
data4 = {{23.5, 6.3195}, {25.5, 6.9645}, {27.5, 7.6667}, {29.5, 
    8.4213}, {31.5, 9.2314}, {33.5, 10.0967}, {35.5, 11.0172}, {37.5, 
    11.9929}, {39.5, 13.0245}, {41.5, 14.11}};

Algorithm

In[13]:= \[Gamma] = 17;
\[Beta] = 0.38;
j1 = 13/2;
j2 = 15/2;
j3 = 4.5;
j4 = 5.5;
j5 = 3.5;
g1 = 2*1/(1 + (5/(16*\[Pi]))^(1/
    2)) (1 - (5/(4*\[Pi]))^(1/2)*\[Beta]*
      Cos[\[Gamma]*\[Pi]/180 + 2/3*\[Pi]*1]);
g2 = 2*1/(1 + (5/(16*\[Pi]))^(1/
    2)) (1 - (5/(4*\[Pi]))^(1/2)*\[Beta]*
      Cos[\[Gamma]*\[Pi]/180 + 2/3*\[Pi]*2]);
g3 = 2*1/(1 + (5/(16*\[Pi]))^(1/
    2)) (1 - (5/(4*\[Pi]))^(1/2)*\[Beta]*
      Cos[\[Gamma]*\[Pi]/180 + 2/3*\[Pi]*3]);
om1[x_, y_] := 
  a*Sqrt[((2 x - 1) (1/g3 - 1/g1) + y/g1)*((2 x - 1) (1/g2 - 1/g1) + y/g1)];
om2[x_, y_] := 
  a*\[Sqrt](((2 y - 1) (1/g3 - 1/g1) + x/g1 + 
        b*(2 y - 1)/(y + 1) Sqrt[
         3] (Sqrt[3] Cos[\[Gamma]*\[Pi]/180] + 
           Sin[\[Gamma]*\[Pi]/180])) ((2 y - 1) (1/g2 - 1/g1) + x/g1 + 
        b*(2 y - 1)/(y + 1) 2 Sqrt[3] Sin[\[Gamma]*\[Pi]/180]));

In[25]:= k1[x_, y_] := (((2 x - 1) (1/g2 - 1/g1) + y/g1)/((2 x - 1) (1/g3 - 1/g1) + y/
     g1)*x^2)^(1/4);
k2[x_, y_] := ((((2 y - 1) (1/g2 - 1/g2) + x/g1 + 
        b*(2 y - 1)/(y (y + 1)) 2 Sqrt[3]
          Sin[\[Gamma]*\[Pi]/180])/((2 y - 1) (1/g3 - 1/g1) + x/g1 + 
        b*(2 y - 1)/(y (y + 1)) Sqrt[
         3] (Sqrt[3] Cos[\[Gamma]*\[Pi]/180] + Sin[\[Gamma]*\[Pi]/180])))*
    y^2)^(1/4);

In[27]:= c1[x_, y_] := 
  om1[x, y]^2 om2[x, y]^2 - 
   a^2*(1/g3^2*k1[x, y]^2*k2[x, y]^2 + 
      x^2*y^2*1/g2^2*1/(k1[x, y]^2*k2[x, y]^2))*om1[x, y]*om2[x, y] + 
   a^4*1/(g2^2*g3^2)*x^2*y^2;
b1[x_, y_] := om1[x, y]^2 + om2[x, y]^2 + 2 a^2*1/(g2*g3)*x*y;

In[29]:= equ[t_, x_, y_] := t^4 - b1[x, y] t^2 + c1[x, y];


In[30]:= p[x_, y_] := Values@Solve[equ[t, x, y] == 0, t];

In[31]:= x1[x_, y_] := p[x, y][[3, 1]];
x2[x_, y_] := p[x, y][[4, 1]];

In[33]:= hmin[x_, y_] := 
  a*((x + y)/2 (1/g2 + 1/g3) + (x^2 + y^2 - x*y)/g1 + 2 x*y*1/g3 - 
     b*(2 y - 1)/(y + 1) Sin[\[Gamma]*\[Pi]/180 + \[Pi]/6]);
h2[x_, y_, n_, m_] := x1[x, y]*(n + 1/2) + x2[x, y]*(m + 1/2);

In[35]:= en1[x_] := hmin[x, j1] + h2[x, j1, 0, 0];
en2[x_] := hmin[x, j1] + h2[x, j1, 1, 0];
en3[x_] := hmin[x, j1] + h2[x, j1, 2, 0];
en4[x_] := hmin[x, j4] + h2[x, j4, 3, 0];

In[39]:= chi2[a_, b_] := 
  Sqrt[((en1[6.5] - tsd1[[1]])^2/tsd1[[1]] + (en1[8.5] - tsd1[[2]])^2/
      tsd1[[2]] + (en1[10.5] - tsd1[[3]])^2/
      tsd1[[3]] + (en1[12.5] - tsd1[[4]])^2/
      tsd1[[4]] + (en1[14.5] - tsd1[[5]])^2/
      tsd1[[5]] + (en1[16.5] - tsd1[[6]])^2/
      tsd1[[6]] + (en1[18.5] - tsd1[[7]])^2/
      tsd1[[7]] + (en1[20.5] - tsd1[[8]])^2/
      tsd1[[8]] + (en1[22.5] - tsd1[[9]])^2/
      tsd1[[9]] + (en1[24.5] - tsd1[[10]])^2/
      tsd1[[10]] + (en1[26.5] - tsd1[[11]])^2/
      tsd1[[11]] + (en1[28.5] - tsd1[[12]])^2/
      tsd1[[12]] + (en1[30.5] - tsd1[[13]])^2/
      tsd1[[13]] + (en1[32.5] - tsd1[[14]])^2/
      tsd1[[14]] + (en1[34.5] - tsd1[[15]])^2/
      tsd1[[15]] + (en1[36.5] - tsd1[[16]])^2/
      tsd1[[16]] + (en1[38.5] - tsd1[[17]])^2/
      tsd1[[17]] + (en1[40.5] - tsd1[[18]])^2/
      tsd1[[18]] + (en1[42.5] - tsd1[[19]])^2/
      tsd1[[19]] + (en1[44.5] - tsd1[[20]])^2/
      tsd1[[20]] + (en1[46.5] - tsd1[[21]])^2/
      tsd1[[21]] + (en1[48.5] - tsd1[[22]])^2/
      tsd1[[22]] + (en2[13.5] - tsd2[[1]])^2/
      tsd2[[1]] + (en2[15.5] - tsd2[[2]])^2/
      tsd2[[2]] + (en2[17.5] - tsd2[[3]])^2/
      tsd2[[3]] + (en2[19.5] - tsd2[[4]])^2/
      tsd2[[4]] + (en2[21.5] - tsd2[[5]])^2/
      tsd2[[5]] + (en2[23.5] - tsd2[[6]])^2/
      tsd2[[6]] + (en2[25.5] - tsd2[[7]])^2/
      tsd2[[7]] + (en2[27.5] - tsd2[[8]])^2/
      tsd2[[8]] + (en2[29.5] - tsd2[[9]])^2/
      tsd2[[9]] + (en2[31.5] - tsd2[[10]])^2/
      tsd2[[10]] + (en2[33.5] - tsd2[[11]])^2/
      tsd2[[11]] + (en2[35.5] - tsd2[[12]])^2/
      tsd2[[12]] + (en2[37.5] - tsd2[[13]])^2/
      tsd2[[13]] + (en2[39.5] - tsd2[[14]])^2/
      tsd2[[14]] + (en2[41.5] - tsd2[[15]])^2/
      tsd2[[15]] + (en2[43.5] - tsd2[[16]])^2/
      tsd2[[16]] + (en2[45.5] - tsd2[[17]])^2/
      tsd2[[17]] + (en3[16.5] - tsd3[[1]])^2/
      tsd3[[1]] + (en3[18.5] - tsd3[[2]])^2/
      tsd3[[2]] + (en3[20.5] - tsd3[[3]])^2/
      tsd3[[3]] + (en3[22.5] - tsd3[[4]])^2/
      tsd3[[4]] + (en3[24.5] - tsd3[[5]])^2/
      tsd3[[5]] + (en3[26.5] - tsd3[[6]])^2/
      tsd3[[6]] + (en3[28.5] - tsd3[[7]])^2/
      tsd3[[7]] + (en3[30.5] - tsd3[[8]])^2/
      tsd3[[8]] + (en3[32.5] - tsd3[[9]])^2/
      tsd3[[10]] + (en3[34.5] - tsd3[[10]])^2/
      tsd3[[10]] + (en3[36.5] - tsd3[[11]])^2/
      tsd3[[11]] + (en3[38.5] - tsd3[[12]])^2/
      tsd3[[12]] + (en3[40.5] - tsd3[[13]])^2/
      tsd3[[13]] + (en3[42.5] - tsd3[[14]])^2/
      tsd3[[14]] + (en4[23.5] - tsd4[[1]])^2/
      tsd4[[1]] + (en4[25.5] - tsd4[[2]])^2/
      tsd4[[2]] + (en4[27.5] - tsd4[[3]])^2/
      tsd4[[3]] + (en4[29.5] - tsd4[[4]])^2/
      tsd4[[4]] + (en4[31.5] - tsd4[[5]])^2/
      tsd4[[5]] + (en4[33.5] - tsd4[[6]])^2/
      tsd4[[6]] + (en4[35.5] - tsd4[[7]])^2/
      tsd4[[7]] + (en4[37.5] - tsd4[[8]])^2/
      tsd4[[8]] + (en4[39.5] - tsd4[[9]])^2/
      tsd4[[9]] + (en4[41.5] - tsd4[[10]])^2/tsd4[[10]])*1/63];
Minimize[chi2[a, b], {a, b}]

Out[40]= If[30477212276628375514 === $SessionID, 
Out[40], Message[
MessageName[Syntax, "noinfoker"]]; Missing["NotAvailable"]; Null]
Attachments:
POSTED BY: Robert Poenaru
6 Replies

So, let's say that I have two functions which depends on two variables, say x and y and the function also has two real parameters a and b. Then you also have a term Q which is a linear combination of the previous functions, and this expression still depends on x and y, with the parameters a and b.

Now you have a final term which is some constant plus the term Q, and this will depend only on x, with a fixed known y.

How would you express such a thing in mathema

POSTED BY: Robert Poenaru

Your description is kind of 'general', i'm not sure what the issue might be. What is the difference between a variable and a parameter? You can use:

f[a_,b_,x_,y_] := …
POSTED BY: Sander Huisman

The expression for:

chi2[a_, b_]

Does not contain a or b…

Simplified example:

f = a^2;
xy[a_] := f
xy[2]

will return a^2, not 4. Because you are using SetDelayed := not Set =

POSTED BY: Sander Huisman

Hello Sander, Thank you for the reply. I modified the := to =. I can get numerical results in chi2.

In[49]:= chi2[1, 2]

Out[49]= 277.249 - 7.44202*10^-16 I

Why should be a problem if in chi2 there is no a and b? The terms a and b are "encoded" in the terms "en[x]" from its expression.

POSTED BY: Robert Poenaru

SetDelayed first 'fills in' a and b with the values you give them in the function, and then evaluates the rhs expression. That is the characteristic of :=, you want first evaluate, then fill in, that is done using =

POSTED BY: Sander Huisman

Compare:

ClearAll[f,xy,a];
f=a^2;
xy[a_]:=f;
xy[2]

ClearAll[f,xy,a];
f=a^2;
xy[a_]=f;
xy[2]

f has a 'encoded' in it. The first and second function definition give different results…

POSTED BY: Sander Huisman
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract