help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

## Re: Leasqr matrix singular to machine precision

 From: Archambault Fabien Subject: Re: Leasqr matrix singular to machine precision Date: Wed, 20 Feb 2008 18:37:20 +0100 User-agent: Thunderbird 2.0.0.4 (X11/20070620)

```Archambault Fabien a écrit :
```
```Hi,

```
I am working on octave since a few days because I need to fit a potential on Lennard-Jones parameters (for those who understand).
```
So I created octave files that can be run in bash.

I think my equations are well written but if possible someone can check.

The script is :

#!/usr/bin/octave -qf
source "octave.script"
x = [ 1.367; 0.224500; 1.768200; -0.120; -0.046000; -0.152100 ];
r(1,1) = RTEMP(1,14);
r(1,2) = RTEMP(2,14);
r(1,3) = RTEMP(3,14);
y(1) = DELTA(14);
r(2,1) = RTEMP(1,1);
r(2,2) = RTEMP(2,1);
r(2,3) = RTEMP(3,1);
y(2) = DELTA(1);
r(3,1) = RTEMP(1,5);
r(3,2) = RTEMP(2,5);
r(3,3) = RTEMP(3,5);
y(3) = DELTA(5);
r(4,1) = RTEMP(1,13);
r(4,2) = RTEMP(2,13);
r(4,3) = RTEMP(3,13);
y(4) = DELTA(13);
r(5,1) = RTEMP(1,9);
r(5,2) = RTEMP(2,9);
r(5,3) = RTEMP(3,9);
y(5) = DELTA(9);
r(6,1) = RTEMP(1,12);
r(6,2) = RTEMP(2,12);
r(6,3) = RTEMP(3,12);
y(6) = DELTA(12);
function y = f(r,x)
```
y = sqrt(x(4)*x(5)) * ( ((x(1)/2+x(2)/2) ./ r(:,1)).^12 - 2 * ((x(1)/2+x(2)/2) ./ r(:,1)).^6 ) \ + sqrt(x(4)*x(5)) * ( ((x(1)/2+x(2)/2) ./ r(:,2)).^12 - 2 * ((x(1)/2+x(2)/2) ./ r(:,2)).^6 ) \ + sqrt(x(4)*x(6)) * ( ((x(1)/2+x(3)/2) ./ r(:,3)).^12 - 2 * ((x(1)/2+x(3)/2) ./ r(:,3)).^6 );
```endfunction
[f1,result,kvg,iter] = leasqr(r,y,x,"f")

See in attach file for octave.script.

When I launch the script the message is :

```
warning: in /usr/share/octave/2.1.73/site/m/octave-forge/optim/leasqr.m near line 299, column 5:
```warning: division by zero
warning: division by zero
```
warning: inverse: matrix singular to machine precision, rcond = 6.93366e-19
```f1 =

-4.8789e-04
-7.2189e-02
-2.5896e-02
-9.8740e-04
-1.0229e-02
-2.1931e-03

result =

1.367000
0.224500
1.768200
-0.120000
-0.046000
-0.152100

kvg = 1
iter = 1

```
The version of octave is 2.1.73 on Mandriva 2007.1. The script leasqr is taken from octave-forge (in rpm files and tested with script available online).
```
```
I have searched also on the mailing list but did not find any clue for this problem.
```
Thanks in advance for your answer,
Fabien Archambault

------------------------------------------------------------------------

_______________________________________________
Help-octave mailing list
address@hidden
https://www.cae.wisc.edu/mailman/listinfo/help-octave
```
```Hi again,

```
I have worked on my scripts since this morning and found some mistakes on my distances so it was impossible to fit with the function.
```
So I modified it :
#!/usr/bin/octave -qf
source "octave.script"
x = [ 1.367; 0.224500; 1.768200; -0.120; -0.046000; -0.152100 ];
r(1,1) = RTEMP(1,8);
r(1,2) = RTEMP(2,8);
r(1,3) = RTEMP(3,8);
y(1) = DELTA(8);
r(2,1) = RTEMP(1,7);
r(2,2) = RTEMP(2,7);
r(2,3) = RTEMP(3,7);
y(2) = DELTA(7);
r(3,1) = RTEMP(1,13);
r(3,2) = RTEMP(2,13);
r(3,3) = RTEMP(3,13);
y(3) = DELTA(13);
r(4,1) = RTEMP(1,2);
r(4,2) = RTEMP(2,2);
r(4,3) = RTEMP(3,2);
y(4) = DELTA(2);
r(5,1) = RTEMP(1,12);
r(5,2) = RTEMP(2,12);
r(5,3) = RTEMP(3,12);
y(5) = DELTA(12);
r(6,1) = RTEMP(1,5);
r(6,2) = RTEMP(2,5);
r(6,3) = RTEMP(3,5);
y(6) = DELTA(5);
function y = f(r,x)
y = \
```
+ sqrt( x(4) * x(5) )*( ( (x(1)+x(2)) ./ r(:,1) ).^12 - 2 * ( (x(1)+x(2)) ./ r(:,1) ).^6) \ + sqrt( x(4) * x(5) )*( ( (x(1)+x(2)) ./ r(:,2) ).^12 - 2 * ( (x(1)+x(2)) ./ r(:,2) ).^6) \ + sqrt( x(4) * x(6) )*( ( (x(1)+x(3)) ./ r(:,3) ).^12 - 2 * ( (x(1)+x(3)) ./ r(:,3) ).^6) \
```;
endfunction
[f1,result,kvg,iter] = leasqr(r,y,x,"f")

And in attach file the new distances.

I have the same error :
```
warning: in /usr/share/octave/2.1.73/site/m/octave-forge/optim/leasqr.m near line 299, column 5:
```warning: division by zero
warning: division by zero
warning: inverse: matrix singular to machine precision, rcond = 6.1851e-23
f1 =

1.66389 +  1.61324i
3.05701 +  1.87630i
-0.06294 +  0.22328i
60.99300 +  2.95108i
-0.10449 +  0.42837i
9.93354 +  2.48499i

result =

1.43557 +   0.00269i
-0.36752 -   1.95093i
1.84242 -   0.00268i
-0.13855 +   0.00670i
204.68952 +   9.96074i
-0.23220 +   0.00222i

kvg = 1
iter = 9

Now it makes iterations but the solutions provided are imaginary !

Do I have to add options like the derivate of the function.

Thanks again for any answer.

--
Fabien Archambault
Equipe de dynamique des assemblages membranaires
Unité Mixte de Recherches CNRS UHP 7565
Université Henri-Poincaré, Nancy I BP 239,
54506 Vandoeuvre-lès-Nancy, cedex France

Tél : 03.83.68.43.96

```
```GROUPANOM1 = [ 20 ];
GROUPAX1 = [ 0 ];
GROUPAY1 = [ 0 ];
GROUPAZ1 = [ 0 ];
GROUPANOM2 = [ 20 ];
GROUPAX2 = [ 0 ];
GROUPAY2 = [ 0 ];
GROUPAZ2 = [ 0 ];
GROUPANOM3 = [ 20 ];
GROUPAX3 = [ 0 ];
GROUPAY3 = [ 0 ];
GROUPAZ3 = [ 0 ];
GROUPANOM4 = [ 20 ];
GROUPAX4 = [ 0 ];
GROUPAY4 = [ 0 ];
GROUPAZ4 = [ 0 ];
GROUPANOM5 = [ 20 ];
GROUPAX5 = [ 0 ];
GROUPAY5 = [ 0 ];
GROUPAZ5 = [ 0 ];
GROUPANOM6 = [ 20 ];
GROUPAX6 = [ 0 ];
GROUPAY6 = [ 0 ];
GROUPAZ6 = [ 0 ];
GROUPANOM7 = [ 20 ];
GROUPAX7 = [ 0 ];
GROUPAY7 = [ 0 ];
GROUPAZ7 = [ 0 ];
GROUPANOM8 = [ 20 ];
GROUPAX8 = [ 0 ];
GROUPAY8 = [ 0 ];
GROUPAZ8 = [ 0 ];
GROUPANOM9 = [ 20 ];
GROUPAX9 = [ 0 ];
GROUPAY9 = [ 0 ];
GROUPAZ9 = [ 0 ];
GROUPANOM10 = [ 20 ];
GROUPAX10 = [ 0 ];
GROUPAY10 = [ 0 ];
GROUPAZ10 = [ 0 ];
GROUPANOM11 = [ 20 ];
GROUPAX11 = [ 0 ];
GROUPAY11 = [ 0 ];
GROUPAZ11 = [ 0 ];
GROUPANOM12 = [ 20 ];
GROUPAX12 = [ 0 ];
GROUPAY12 = [ 0 ];
GROUPAZ12 = [ 0 ];
GROUPANOM13 = [ 20 ];
GROUPAX13 = [ 0 ];
GROUPAY13 = [ 0 ];
GROUPAZ13 = [ 0 ];
GROUPANOM14 = [ 20 ];
GROUPAX14 = [ 0 ];
GROUPAY14 = [ 0 ];
GROUPAZ14 = [ 0 ];
GROUPANOM15 = [ 20 ];
GROUPAX15 = [ 0 ];
GROUPAY15 = [ 0 ];
GROUPAZ15 = [ 0 ];
GROUPBNOM1 = [ 1,1,8 ];
GROUPBX1 = [ -.756358,.756358,0 ];
GROUPBY1 = [ 0,0,0 ];
GROUPBZ1 = [ -2.588115,-2.588115,-1.900000 ];
GROUPBNOM2 = [ 1,1,8 ];
GROUPBX2 = [ -.756358,.756358,0 ];
GROUPBY2 = [ 0,0,0 ];
GROUPBZ2 = [
-2.68811500000000000000,-2.68811500000000000000,-2.00000000000000000000 ];
GROUPBNOM3 = [ 1,1,8 ];
GROUPBX3 = [ -.756358,.756358,0 ];
GROUPBY3 = [ 0,0,0 ];
GROUPBZ3 = [
-2.78811500000000000000,-2.78811500000000000000,-2.10000000000000000000 ];
GROUPBNOM4 = [ 1,1,8 ];
GROUPBX4 = [ -.756358,.756358,0 ];
GROUPBY4 = [ 0,0,0 ];
GROUPBZ4 = [
-2.88811500000000000000,-2.88811500000000000000,-2.20000000000000000000 ];
GROUPBNOM5 = [ 1,1,8 ];
GROUPBX5 = [ -.756358,.756358,0 ];
GROUPBY5 = [ 0,0,0 ];
GROUPBZ5 = [
-2.98811500000000000000,-2.98811500000000000000,-2.30000000000000000000 ];
GROUPBNOM6 = [ 1,1,8 ];
GROUPBX6 = [ -.756358,.756358,0 ];
GROUPBY6 = [ 0,0,0 ];
GROUPBZ6 = [
-3.08811500000000000000,-3.08811500000000000000,-2.40000000000000000000 ];
GROUPBNOM7 = [ 1,1,8 ];
GROUPBX7 = [ -.756358,.756358,0 ];
GROUPBY7 = [ 0,0,0 ];
GROUPBZ7 = [
-3.18811500000000000000,-3.18811500000000000000,-2.50000000000000000000 ];
GROUPBNOM8 = [ 1,1,8 ];
GROUPBX8 = [ -.756358,.756358,0 ];
GROUPBY8 = [ 0,0,0 ];
GROUPBZ8 = [
-3.28811500000000000000,-3.28811500000000000000,-2.60000000000000000000 ];
GROUPBNOM9 = [ 1,1,8 ];
GROUPBX9 = [ -.756358,.756358,0 ];
GROUPBY9 = [ 0,0,0 ];
GROUPBZ9 = [
-3.38811500000000000000,-3.38811500000000000000,-2.70000000000000000000 ];
GROUPBNOM10 = [ 1,1,8 ];
GROUPBX10 = [ -.756358,.756358,0 ];
GROUPBY10 = [ 0,0,0 ];
GROUPBZ10 = [
-3.48811500000000000000,-3.48811500000000000000,-2.80000000000000000000 ];
GROUPBNOM11 = [ 1,1,8 ];
GROUPBX11 = [ -.756358,.756358,0 ];
GROUPBY11 = [ 0,0,0 ];
GROUPBZ11 = [
-3.58811500000000000000,-3.58811500000000000000,-2.90000000000000000000 ];
GROUPBNOM12 = [ 1,1,8 ];
GROUPBX12 = [ -.756358,.756358,0 ];
GROUPBY12 = [ 0,0,0 ];
GROUPBZ12 = [
-4.18811500000000000000,-4.18811500000000000000,-3.50000000000000000000 ];
GROUPBNOM13 = [ 1,1,8 ];
GROUPBX13 = [ -.756358,.756358,0 ];
GROUPBY13 = [ 0,0,0 ];
GROUPBZ13 = [
-4.68811500000000000000,-4.68811500000000000000,-4.00000000000000000000 ];
GROUPBNOM14 = [ 1,1,8 ];
GROUPBX14 = [ -.756358,.756358,0 ];
GROUPBY14 = [ 0,0,0 ];
GROUPBZ14 = [
-5.18811500000000000000,-5.18811500000000000000,-4.50000000000000000000 ];
GROUPBNOM15 = [ 1,1,8 ];
GROUPBX15 = [ -.756358,.756358,0 ];
GROUPBY15 = [ 0,0,0 ];
GROUPBZ15 = [
-5.68811500000000000000,-5.68811500000000000000,-5.00000000000000000000 ];
DELTA = [
100.03,65.9,43.63,28.9,19.07,12.48,8.05,5.07,3.08,1.76,0.9,-0.34,-0.2,-0.09,-0.04
];
RTEMP =
[sqrt((GROUPAX1(1)-1*GROUPBX1(1))^2+(GROUPAY1(1)-1*GROUPBY1(1))^2+(GROUPAZ1(1)-1*GROUPBZ1(1))^2),sqrt((GROUPAX2(1)-1*GROUPBX2(1))^2+(GROUPAY2(1)-1*GROUPBY2(1))^2+(GROUPAZ2(1)-1*GROUPBZ2(1))^2),sqrt((GROUPAX3(1)-1*GROUPBX3(1))^2+(GROUPAY3(1)-1*GROUPBY3(1))^2+(GROUPAZ3(1)-1*GROUPBZ3(1))^2),sqrt((GROUPAX4(1)-1*GROUPBX4(1))^2+(GROUPAY4(1)-1*GROUPBY4(1))^2+(GROUPAZ4(1)-1*GROUPBZ4(1))^2),sqrt((GROUPAX5(1)-1*GROUPBX5(1))^2+(GROUPAY5(1)-1*GROUPBY5(1))^2+(GROUPAZ5(1)-1*GROUPBZ5(1))^2),sqrt((GROUPAX6(1)-1*GROUPBX6(1))^2+(GROUPAY6(1)-1*GROUPBY6(1))^2+(GROUPAZ6(1)-1*GROUPBZ6(1))^2),sqrt((GROUPAX7(1)-1*GROUPBX7(1))^2+(GROUPAY7(1)-1*GROUPBY7(1))^2+(GROUPAZ7(1)-1*GROUPBZ7(1))^2),sqrt((GROUPAX8(1)-1*GROUPBX8(1))^2+(GROUPAY8(1)-1*GROUPBY8(1))^2+(GROUPAZ8(1)-1*GROUPBZ8(1))^2),sqrt((GROUPAX9(1)-1*GROUPBX9(1))^2+(GROUPAY9(1)-1*GROUPBY9(1))^2+(GROUPAZ9(1)-1*GROUPBZ9(1))^2),sqrt((GROUPAX10(1)-1*GROUPBX10(1))^2+(GROUPAY10(1)-1*GROUPBY10(1))^2+(GROUPAZ10(1)-1*GROUPBZ10(1))^2),sqrt((GROUPAX11(1)-1*GROUPBX11(1))^2+(GROUPAY11(1)-1*GROUPBY11(1))^2+(GROUPAZ11(1)-1*GROUPBZ11(1))^2),sqrt((GROUPAX12(1)-1*GROUPBX12(1))^2+(GROUPAY12(1)-1*GROUPBY12(1))^2+(GROUPAZ12(1)-1*GROUPBZ12(1))^2),sqrt((GROUPAX13(1)-1*GROUPBX13(1))^2+(GROUPAY13(1)-1*GROUPBY13(1))^2+(GROUPAZ13(1)-1*GROUPBZ13(1))^2),sqrt((GROUPAX14(1)-1*GROUPBX14(1))^2+(GROUPAY14(1)-1*GROUPBY14(1))^2+(GROUPAZ14(1)-1*GROUPBZ14(1))^2),sqrt((GROUPAX15(1)-1*GROUPBX15(1))^2+(GROUPAY15(1)-1*GROUPBY15(1))^2+(GROUPAZ15(1)-1*GROUPBZ15(1))^2);sqrt((GROUPAX1(1)-1*GROUPBX1(2))^2+(GROUPAY1(1)-1*GROUPBY1(2))^2+(GROUPAZ1(1)-1*GROUPBZ1(2))^2),sqrt((GROUPAX2(1)-1*GROUPBX2(2))^2+(GROUPAY2(1)-1*GROUPBY2(2))^2+(GROUPAZ2(1)-1*GROUPBZ2(2))^2),sqrt((GROUPAX3(1)-1*GROUPBX3(2))^2+(GROUPAY3(1)-1*GROUPBY3(2))^2+(GROUPAZ3(1)-1*GROUPBZ3(2))^2),sqrt((GROUPAX4(1)-1*GROUPBX4(2))^2+(GROUPAY4(1)-1*GROUPBY4(2))^2+(GROUPAZ4(1)-1*GROUPBZ4(2))^2),sqrt((GROUPAX5(1)-1*GROUPBX5(2))^2+(GROUPAY5(1)-1*GROUPBY5(2))^2+(GROUPAZ5(1)-1*GROUPBZ5(2))^2),sqrt((GROUPAX6(1)-1*GROUPBX6(2))^2+(GROUPAY6(1)-1*GROUPBY6(2))^2+(GROUPAZ6(1)-1*GROUPBZ6(2))^2),sqrt((GROUPAX7(1)-1*GROUPBX7(2))^2+(GROUPAY7(1)-1*GROUPBY7(2))^2+(GROUPAZ7(1)-1*GROUPBZ7(2))^2),sqrt((GROUPAX8(1)-1*GROUPBX8(2))^2+(GROUPAY8(1)-1*GROUPBY8(2))^2+(GROUPAZ8(1)-1*GROUPBZ8(2))^2),sqrt((GROUPAX9(1)-1*GROUPBX9(2))^2+(GROUPAY9(1)-1*GROUPBY9(2))^2+(GROUPAZ9(1)-1*GROUPBZ9(2))^2),sqrt((GROUPAX10(1)-1*GROUPBX10(2))^2+(GROUPAY10(1)-1*GROUPBY10(2))^2+(GROUPAZ10(1)-1*GROUPBZ10(2))^2),sqrt((GROUPAX11(1)-1*GROUPBX11(2))^2+(GROUPAY11(1)-1*GROUPBY11(2))^2+(GROUPAZ11(1)-1*GROUPBZ11(2))^2),sqrt((GROUPAX12(1)-1*GROUPBX12(2))^2+(GROUPAY12(1)-1*GROUPBY12(2))^2+(GROUPAZ12(1)-1*GROUPBZ12(2))^2),sqrt((GROUPAX13(1)-1*GROUPBX13(2))^2+(GROUPAY13(1)-1*GROUPBY13(2))^2+(GROUPAZ13(1)-1*GROUPBZ13(2))^2),sqrt((GROUPAX14(1)-1*GROUPBX14(2))^2+(GROUPAY14(1)-1*GROUPBY14(2))^2+(GROUPAZ14(1)-1*GROUPBZ14(2))^2),sqrt((GROUPAX15(1)-1*GROUPBX15(2))^2+(GROUPAY15(1)-1*GROUPBY15(2))^2+(GROUPAZ15(1)-1*GROUPBZ15(2))^2);sqrt((GROUPAX1(1)-1*GROUPBX1(3))^2+(GROUPAY1(1)-1*GROUPBY1(3))^2+(GROUPAZ1(1)-1*GROUPBZ1(3))^2),sqrt((GROUPAX2(1)-1*GROUPBX2(3))^2+(GROUPAY2(1)-1*GROUPBY2(3))^2+(GROUPAZ2(1)-1*GROUPBZ2(3))^2),sqrt((GROUPAX3(1)-1*GROUPBX3(3))^2+(GROUPAY3(1)-1*GROUPBY3(3))^2+(GROUPAZ3(1)-1*GROUPBZ3(3))^2),sqrt((GROUPAX4(1)-1*GROUPBX4(3))^2+(GROUPAY4(1)-1*GROUPBY4(3))^2+(GROUPAZ4(1)-1*GROUPBZ4(3))^2),sqrt((GROUPAX5(1)-1*GROUPBX5(3))^2+(GROUPAY5(1)-1*GROUPBY5(3))^2+(GROUPAZ5(1)-1*GROUPBZ5(3))^2),sqrt((GROUPAX6(1)-1*GROUPBX6(3))^2+(GROUPAY6(1)-1*GROUPBY6(3))^2+(GROUPAZ6(1)-1*GROUPBZ6(3))^2),sqrt((GROUPAX7(1)-1*GROUPBX7(3))^2+(GROUPAY7(1)-1*GROUPBY7(3))^2+(GROUPAZ7(1)-1*GROUPBZ7(3))^2),sqrt((GROUPAX8(1)-1*GROUPBX8(3))^2+(GROUPAY8(1)-1*GROUPBY8(3))^2+(GROUPAZ8(1)-1*GROUPBZ8(3))^2),sqrt((GROUPAX9(1)-1*GROUPBX9(3))^2+(GROUPAY9(1)-1*GROUPBY9(3))^2+(GROUPAZ9(1)-1*GROUPBZ9(3))^2),sqrt((GROUPAX10(1)-1*GROUPBX10(3))^2+(GROUPAY10(1)-1*GROUPBY10(3))^2+(GROUPAZ10(1)-1*GROUPBZ10(3))^2),sqrt((GROUPAX11(1)-1*GROUPBX11(3))^2+(GROUPAY11(1)-1*GROUPBY11(3))^2+(GROUPAZ11(1)-1*GROUPBZ11(3))^2),sqrt((GROUPAX12(1)-1*GROUPBX12(3))^2+(GROUPAY12(1)-1*GROUPBY12(3))^2+(GROUPAZ12(1)-1*GROUPBZ12(3))^2),sqrt((GROUPAX13(1)-1*GROUPBX13(3))^2+(GROUPAY13(1)-1*GROUPBY13(3))^2+(GROUPAZ13(1)-1*GROUPBZ13(3))^2),sqrt((GROUPAX14(1)-1*GROUPBX14(3))^2+(GROUPAY14(1)-1*GROUPBY14(3))^2+(GROUPAZ14(1)-1*GROUPBZ14(3))^2),sqrt((GROUPAX15(1)-1*GROUPBX15(3))^2+(GROUPAY15(1)-1*GROUPBY15(3))^2+(GROUPAZ15(1)-1*GROUPBZ15(3))^2)
];
```

reply via email to

 [Prev in Thread] Current Thread [Next in Thread]