| Author |
 Topic  |
|
|
DrKlein1
Germany
6 Posts |
 Posted - 11/23/2010 : 07:45:35 AM
|
Origin Ver. 6.1 v6.1052(B232) and Service Release (Select Help-->About Origin):
Operating System: Windows XP
I would like to use the following script to iteratively calculate n and k from a starting guess. The calculation is not complicated but nested.
I think due to a lack of precision, the returned values from Origin do not converge:
n=3.5; //starting guess
k=0.1; //starting guess
d=80;
ns=1.45;
l=1900;
R=1.29744976541021;
R'=1.29475398678337;
loop(i,1,5){
r=(1+ns^2)*(n^2+k^2) - (n^2+k^2)^2 -ns^2-4*ns*k^2;
S=2*k*(ns-1)*(n^2+k^2+ns);
P=2*( (n^2+k^2+ns^2)*(n^2+k^2+1) - 4*n^2*ns);
Q=4*n*( (n^2+k^2+1)*ns - (n^2+k^2+ns^2) );
alpha=2*pi*k*d/l;
gamma=2*pi*n*d/l;
Pn=4*n*((n^2 + k^2 + 1) + (n^2 + k^2 + ns^2) - 4*ns);
Pk=4*k*((n^2 + k^2 + 1) + (n^2 + k^2 + ns^2));
Qn=4*((n^2 + k^2 + 1)*ns - (n^2 + k^2 + ns^2)) + 8*n^2*(ns-1);
Qk=8*n*k*(ns-1);
Rn=-2*n*((n^2 + k^2 - 1) + (n^2 + k^2 - ns^2));
Rk=-2*k*((n^2 + k^2 - 1) + (n^2 + k^2 - ns^2) + 4*ns);
Sn=4*n*k*(ns-1);
Sk=2*(ns-1)*(n^2 + k^2 + ns) + 4*k^2*(ns-1);
f1=P*cosh(2*alpha)+Q*sinh(2*alpha)+2*r*cos(2*gamma)+2*S*sin(2*gamma)-16*ns*(n^2 + k^2) * (R);
f2=P*cosh(2*alpha)-Q*sinh(2*alpha)+2*r*cos(2*gamma)-2*S*sin(2*gamma)-16*ns*(n^2 + k^2) * (R');
f1n=Pn*cosh(2*alpha)+Qn*sinh(2*alpha)+2*Rn*cos(2*gamma)+2*Sn*sin(2*gamma)+(8*pi*d/l) * (S*cos(2*gamma)-r*sin(2*gamma))-32*ns*n * (R);
f2n=Pn*cosh(2*alpha)-Qn*sinh(2*alpha)+2*Rn*cos(2*gamma)-2*Sn*sin(2*gamma)-(8*pi*d/l) * (S*cos(2*gamma)+r*sin(2*gamma))-32*ns*n * (R');
f1k=Pk*cosh(2*alpha)+Qk*sinh(2*alpha)+2*Rk*cos(2*gamma)+2*Sk*sin(2*gamma)+(4*pi*d/l) * (P*sinh(2*alpha)+Q*cosh(2*alpha))-32*ns*k * (R);
f2k=Pk*cosh(2*alpha)-Qk*sinh(2*alpha)+2*Rk*cos(2*gamma)-2*Sk*sin(2*gamma)+(4*pi*d/l) * (P*sinh(2*alpha)-Q*cosh(2*alpha))-32*ns*k * (R');
n=n + (f2*f1k - f1*f2k) / (f1n*f2k - f1k*f2n); //recursive formula
k=k + (f1*f2n - f2*f1n) / (f1n*f2k - f1k*f2n); //recursive formula
n=;
k=;
}
//this is what Origin returns:
N=2,115882
K=-10,41146
N=-39,81687
K=-13,19948
N=-17,0642
K=-15,14635
N=-0,1351484
K=-17,03149
N=98,40254
K=-15,56333
//this is what Mathematica returns, {n,k}=
3.72931 -0.00867624
3.5835 0.004228
3.64824 0.00609837
3.61006 0.00583065
3.63044 0.00594162
3.61885 0.00587884
3.62523 0.00591373
3.62165 0.00589426
3.62364 0.0059051
3.62253 0.00589906
3.62315 0.00590243
3.62281 0.00590055
3.623 0.00590159
3.62289 0.00590101
3.62295 0.00590134
3.62292 0.00590116
3.62294 0.00590126
3.62293 0.0059012
3.62293 0.00590123
Is there a way to increase precision, so that it works in Origin, too ? |
|
|
larry_lan
China
Posts |
Posted - 11/24/2010 : 03:07:23 AM
|
Hi:
LabTalk numeric variables hold double precision real numbers, which consistant with C/C++ double type. The main problem in your script is, LabTalk variable is case-insensitive, while Mathematica is case-sensitive. Note the variable r & R! They are the same in LabTalk, but different in Mathematica. I think you can get resets which close to Mathematica by the following script (I changed your r into rr, and changed your R' into R2).
n=3.5;
k=0.1;
d=80;
ns=1.45;
l=1900;
R=1.29744976541021;
R2=1.29475398678337;
loop(ii,1,5){
rr=(1+ns^2)*(n^2+k^2) - (n^2+k^2)^2 -ns^2-4*ns*k^2;
S=2*k*(ns-1)*(n^2+k^2+ns);
P=2*( (n^2+k^2+ns^2)*(n^2+k^2+1) - 4*n^2*ns);
Q=4*n*( (n^2+k^2+1)*ns - (n^2+k^2+ns^2) );
alpha=2*pi*k*d/l;
gamma=2*pi*n*d/l;
Pn=4*n*((n^2 + k^2 + 1) + (n^2 + k^2 + ns^2) - 4*ns);
Pk=4*k*((n^2 + k^2 + 1) + (n^2 + k^2 + ns^2));
Qn=4*((n^2 + k^2 + 1)*ns - (n^2 + k^2 + ns^2)) + 8*n^2*(ns-1);
Qk=8*n*k*(ns-1);
Rn=-2*n*((n^2 + k^2 - 1) + (n^2 + k^2 - ns^2));
Rk=-2*k*((n^2 + k^2 - 1) + (n^2 + k^2 - ns^2) + 4*ns);
Sn=4*n*k*(ns-1);
Sk=2*(ns-1)*(n^2 + k^2 + ns) + 4*k^2*(ns-1);
f1=P*cosh(2*alpha)+Q*sinh(2*alpha)+2*rr*cos(2*gamma)+2*S*sin(2*gamma)-16*ns*(n^2 + k^2) * (R);
f2=P*cosh(2*alpha)-Q*sinh(2*alpha)+2*rr*cos(2*gamma)-2*S*sin(2*gamma)-16*ns*(n^2 + k^2) * (R2);
f1n=Pn*cosh(2*alpha)+Qn*sinh(2*alpha)+2*Rn*cos(2*gamma)+2*Sn*sin(2*gamma)+(8*pi*d/l) * (S*cos(2*gamma)-rr*sin(2*gamma))-32*ns*n * (R);
f2n=Pn*cosh(2*alpha)-Qn*sinh(2*alpha)+2*Rn*cos(2*gamma)-2*Sn*sin(2*gamma)-(8*pi*d/l) * (S*cos(2*gamma)+rr*sin(2*gamma))-32*ns*n * (R2);
f1k=Pk*cosh(2*alpha)+Qk*sinh(2*alpha)+2*Rk*cos(2*gamma)+2*Sk*sin(2*gamma)+(4*pi*d/l) * (P*sinh(2*alpha)+Q*cosh(2*alpha))-32*ns*k * (R);
f2k=Pk*cosh(2*alpha)-Qk*sinh(2*alpha)+2*Rk*cos(2*gamma)-2*Sk*sin(2*gamma)+(4*pi*d/l) * (P*sinh(2*alpha)-Q*cosh(2*alpha))-32*ns*k * (R2);
n=n + (f2*f1k - f1*f2k) / (f1n*f2k - f1k*f2n);
k=k + (f1*f2n - f2*f1n) / (f1n*f2k - f1k*f2n);
ty "Round $(ii), n = $(n), k = $(k)";
}
Here is Origin output:
Round 1, n = 3.7293134732918, k = -0.0086762444755236
Round 2, n = 3.6905812477565, k = 0.0058943423583641
Round 3, n = 3.6892535795924, k = 0.0062688087745428
Round 4, n = 3.6892520775642, k = 0.0062691354421516
Round 5, n = 3.6892520775623, k = 0.0062691354424736
And n is converged to 3.68925.
In Mathematica, if you run
n=3.5;
k=0.1;
d=80;
ns=1.45;
l=1900;
R=1.29744976541021;
R2=1.29475398678337;
For[i=1, i<=5, i++,
r=(1+ns^2)*(n^2+k^2) - (n^2+k^2)^2 -ns^2-4*ns*k^2;
S=2*k*(ns-1)*(n^2+k^2+ns);
P=2*( (n^2+k^2+ns^2)*(n^2+k^2+1) - 4*n^2*ns);
Q=4*n*( (n^2+k^2+1)*ns - (n^2+k^2+ns^2) );
alpha=2*\[Pi]*k*d/l;
gamma=2*\[Pi]*n*d/l;
Pn=4*n*((n^2 + k^2 + 1) + (n^2 + k^2 + ns^2) - 4*ns);
Pk=4*k*((n^2 + k^2 + 1) + (n^2 + k^2 + ns^2));
Qn=4*((n^2 + k^2 + 1)*ns - (n^2 + k^2 + ns^2)) + 8*n^2*(ns-1);
Qk=8*n*k*(ns-1);
Rn=-2*n*((n^2 + k^2 - 1) + (n^2 + k^2 - ns^2));
Rk=-2*k*((n^2 + k^2 - 1) + (n^2 + k^2 - ns^2) + 4*ns);
Sn=4*n*k*(ns-1);
Sk=2*(ns-1)*(n^2 + k^2 + ns) + 4*k^2*(ns-1);
f1=P*Cosh[2*alpha]+Q*Sinh[2*alpha]+2*r*Cos[2*gamma]+2*S*Sin[2*gamma]-16*ns*(n^2 + k^2) * (R);
f2=P*Cosh[2*alpha]-Q*Sinh[2*alpha]+2*r*Cos[2*gamma]-2*S*Sin[2*gamma]-16*ns*(n^2 + k^2) * (R2);
f1n=Pn*Cosh[2*alpha]+Qn*Sinh[2*alpha]+2*Rn*Cos[2*gamma]+2*Sn*Sin[2*gamma]+(8*\[Pi]*d/l) * (S*Cos[2*gamma]-r*Sin[2*gamma])-32*ns*n * (R);
f2n=Pn*Cosh[2*alpha]-Qn*Sinh[2*alpha]+2*Rn*Cos[2*gamma]-2*Sn*Sin[2*gamma]-(8*\[Pi]*d/l) * (S*Cos[2*gamma]+r*Sin[2*gamma])-32*ns*n * (R2);
f1k=Pk*Cosh[2*alpha]+Qk*Sinh[2*alpha]+2*Rk*Cos[2*gamma]+2*Sk*Sin[2*gamma]+(4*\[Pi]*d/l) * (P*Sinh[2*alpha]+Q*Cosh[2*alpha])-32*ns*k * (R);
f2k=Pk*Cosh[2*alpha]-Qk*Sinh[2*alpha]+2*Rk*Cos[2*gamma]-2*Sk*Sin[2*gamma]+(4*\[Pi]*d/l) * (P*Sinh[2*alpha]-Q*Cosh[2*alpha])-32*ns*k * (R2);
n=n + (f2*f1k - f1*f2k) / (f1n*f2k - f1k*f2n);
k=k + (f1*f2n - f2*f1n) / (f1n*f2k - f1k*f2n);
Print["n = ", n, " and k = ", k];
]
You will get:
n = 3.72931 and k = -0.00867624
n = 3.69058 and k = 0.00589434
n = 3.68925 and k = 0.00626881
n = 3.68925 and k = 0.00626914
n = 3.68925 and k = 0.00626914
You can see the results are consistent. And maybe you need to set WorkingPrecision or PrecisionGoal in Mathematica to see more precision?
I don't know why your Mathematica output 3.62293, maybe you can debug your Mathematica code by yourself.
Thanks
Larry
OriginLab |
Edited by - larry_lan on 11/24/2010 04:04:37 AM |
 |
|
|
DrKlein1
Germany
6 Posts |
Posted - 11/24/2010 : 09:22:41 AM
|
| Thank you very much ! |
|
|
| |
Topic |
|
|
|
Increase numeric precision