[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)
];
- Leasqr matrix singular to machine precision, Archambault Fabien, 2008/02/20
- Re: Leasqr matrix singular to machine precision,
Archambault Fabien <=
- Re: Leasqr matrix singular to machine precision, Archambault Fabien, 2008/02/21
- Re: Leasqr matrix singular to machine precision, Francesco Potorti`, 2008/02/21
- Re: Leasqr matrix singular to machine precision, Archambault Fabien, 2008/02/22
- Re: Leasqr matrix singular to machine precision, Francesco Potorti`, 2008/02/22
- Re: Leasqr matrix singular to machine precision, Archambault Fabien, 2008/02/22
- Re: Leasqr matrix singular to machine precision, Przemek Klosowski, 2008/02/25
- Re: Leasqr matrix singular to machine precision, Archambault Fabien, 2008/02/25
- Re: Leasqr matrix singular to machine precision, Dmitri A. Sergatskov, 2008/02/25
- Re: Leasqr matrix singular to machine precision, Archambault Fabien, 2008/02/26
- Re: Leasqr matrix singular to machine precision, Doug Stewart, 2008/02/26