octave:2> catenaria1 octave:3> catenaria1 ans = 2.2540 ans = 1.1843 ans = 0.49160 ans = 0.15354 ans = 0.041351 ans = 0.0097118 ans = 0.0020674 ans = 4.0616e-04 ans = 2.7499 ans = 1.4131 ans = 0.56371 ans = 0.16622 ans = 0.041323 ans = 0.0087770 ans = 0.0016538 ans = 2.8175e-04 ans = 4.4126e-05 ans = 3.1710 ans = 1.6170 ans = 0.63599 ans = 0.18368 ans = 0.044336 ans = 0.0090647 ans = 0.0016286 ans = 2.6203e-04 ans = 3.8378e-05 ans = 3.5430 ans = 1.8002 ans = 0.70346 ans = 0.20119 ans = 0.047891 ans = 0.0096144 ans = 0.0016879 ans = 2.6404e-04 ans = 3.7399e-05 ans = 3.8799 ans = 1.9676 ans = 0.76611 ans = 0.21794 ans = 0.051483 ans = 0.010232 ans = 0.0017734 ans = 2.7306e-04 ans = 3.7951e-05 octave:4> catenaria1 ans = 1.4144 ans = 0.027983 ans = 3.6796e-06 ans = 1.7359 ans = 0.033794 ans = 4.1921e-06 ans = 2.0060 ans = 0.038831 ans = 4.7191e-06 ans = 2.2435 ans = 0.043315 ans = 5.2143e-06 ans = 2.4581 ans = 0.047391 ans = 5.6756e-06 octave:5> errori=[1.4144,0.027983,3.6796e-06] errori = 1.4144e+00 2.7983e-02 3.6796e-06 octave:6> errori(2)/errori(1)^2 ans = 0.013988 octave:7> errori(3)/errori(2)^2 ans = 0.0046991 octave:8> F F = @(u) A * u - a * sqrt (1 + (B * u) .^ 2) - [1 / h ^ 2 - a; zeros(m - 2, 1); 1 / h ^ 2 - a] octave:9> J J = @(u) A - a * spdiags (B * u ./ sqrt (1 + (B * u) .^ 2), 0, m, m) * B octave:10> u u = 1.000000 0.883195 0.774237 0.672640 0.577954 0.489756 0.407654 0.331283 0.260303 0.194398 0.133275 0.076662 0.024307 -0.024023 -0.068542 -0.109450 -0.146927 -0.181141 -0.212244 -0.240375 -0.265657 -0.288205 -0.308119 -0.325486 -0.340384 -0.352880 -0.363029 -0.370876 -0.376456 -0.379795 -0.380906 -0.379795 -0.376456 -0.370876 -0.363029 -0.352880 -0.340384 -0.325486 -0.308119 -0.288205 -0.265657 -0.240375 -0.212244 -0.181141 -0.146927 -0.109450 -0.068542 -0.024023 0.024307 0.076662 0.133275 0.194398 0.260303 0.331283 0.407654 0.489756 0.577954 0.672640 0.774237 0.883195 1.000000 octave:11> u=rand(size(u0)); octave:12> v=rand(size(u0)); octave:13> J(u)*v ans = 235.7587 -5.4523 21.4914 -238.1317 860.6682 -1410.1716 1217.8192 -868.1574 752.4471 -83.4672 -50.5633 -1007.1955 1247.6598 -598.9449 224.0251 -521.3904 640.0942 -81.7299 -231.0231 -139.1602 717.2188 -880.1316 1082.1103 -774.4509 -122.6550 24.8397 343.0723 -201.2472 253.9453 13.7738 -512.5255 279.3462 406.7082 -310.5327 -201.0418 -190.0057 506.1570 -524.9669 1121.8627 -1055.7330 -225.5126 572.0958 195.8739 -203.2273 -735.1937 1088.2298 -439.8002 -52.5301 -28.0329 275.2496 -639.0120 1030.7310 -889.6385 -23.8195 481.8739 250.7220 -459.1462 -38.1552 35.5962 -40.6024 797.2820 octave:14> octave:14> octave:14> octave:14> octave:14> octave:14> octave:14> octave:14> octave:14> octave:14> octave:14> octave:14> epsilon=1e-4 epsilon = 1.0000e-04 octave:15> (F(u+epsilon*v)-F(u))/epsilon ans = 235.7587 -5.4523 21.4914 -238.1317 860.6670 -1410.1716 1217.8192 -868.1574 752.4470 -83.4674 -50.5635 -1007.1955 1247.6569 -598.9453 224.0251 -521.3911 640.0942 -81.7299 -231.0231 -139.1602 717.2188 -880.1316 1082.1103 -774.4512 -122.6555 24.8396 343.0723 -201.2472 253.9453 13.7738 -512.5255 279.3462 406.7082 -310.5370 -201.0418 -190.0057 506.1569 -524.9669 1121.8627 -1055.7331 -225.5135 572.0956 195.8736 -203.2273 -735.1939 1088.2298 -439.8002 -52.5301 -28.0329 275.2495 -639.0120 1030.7309 -889.6386 -23.8196 481.8739 250.7220 -459.1511 -38.1552 35.5961 -40.6024 797.2820 octave:16> w1=J(u)*v; octave:17> w2=(F(u+epsilon*v)-F(u))/epsilon; octave:18> norm(w1-w2) ans = 0.0073022 octave:19> catenaria1 ans = 2.2540 ans = 1.1843 ans = 0.49160 ans = 0.15354 ans = 0.041351 ans = 0.0097118 ans = 0.0020674 ans = 4.0616e-04 ans = 2.7499 ans = 1.4131 ans = 0.56371 ans = 0.16622 ans = 0.041323 ans = 0.0087770 ans = 0.0016538 ans = 2.8175e-04 ans = 4.4126e-05 ans = 3.1710 ans = 1.6170 ans = 0.63599 ans = 0.18368 ans = 0.044336 ans = 0.0090647 ans = 0.0016286 ans = 2.6203e-04 ans = 3.8378e-05 ans = 3.5430 ans = 1.8002 ans = 0.70346 ans = 0.20119 ans = 0.047891 ans = 0.0096144 ans = 0.0016879 ans = 2.6404e-04 ans = 3.7399e-05 ans = 3.8799 ans = 1.9676 ans = 0.76611 ans = 0.21794 ans = 0.051483 ans = 0.010232 ans = 0.0017734 ans = 2.7306e-04 ans = 3.7951e-05 octave:20> octave:20> u=rand(size(u0)); octave:21> v=rand(size(u0)); octave:22> w1=J(u)*v; octave:23> w2=(F(u+epsilon*v)-F(u))/epsilon; error: `epsilon' undefined near line 23 column 9 error: evaluating argument list element number 1 octave:23> epsilon=1e-4 epsilon = 1.0000e-04 octave:24> w2=(F(u+epsilon*v)-F(u))/epsilon; octave:25> norm(w1-w2) ans = 93.544 octave:26> catenaria1 ans = 0.96817 ans = 0.34334 ans = 0.10042 ans = 0.024839 ans = 0.0053565 ans = 0.0010332 ans = 1.8191e-04 ans = 1.1900 ans = 0.42153 ans = 0.12222 ans = 0.029651 ans = 0.0061887 ans = 0.0011380 ans = 1.8796e-04 ans = 1.3758 ans = 0.48716 ans = 0.14081 ans = 0.033924 ans = 0.0069960 ans = 0.0012635 ans = 2.0356e-04 ans = 2.9694e-05 ans = 1.5391 ans = 0.54487 ans = 0.15727 ans = 0.037765 ans = 0.0077442 ans = 0.0013866 ans = 2.2073e-04 ans = 3.1692e-05 ans = 1.6865 ans = 0.59700 ans = 0.17218 ans = 0.041272 ans = 0.0084372 ans = 0.0015036 ans = 2.3775e-04 ans = 3.3835e-05 octave:27> octave:27> mrif=1001 mrif = 1001 octave:28> m=51 m = 51 octave:29> 1:(mrif-1)/(m-1):1001 ans = Columns 1 through 11: 1 21 41 61 81 101 121 141 161 181 201 Columns 12 through 22: 221 241 261 281 301 321 341 361 381 401 421 Columns 23 through 33: 441 461 481 501 521 541 561 581 601 621 641 Columns 34 through 44: 661 681 701 721 741 761 781 801 821 841 861 Columns 45 through 51: 881 901 921 941 961 981 1001 octave:30> mrif=1000 mrif = 1000 octave:31> m=20 m = 20 octave:32> 1:(mrif-1)/(m-1):1000 ans = Columns 1 through 5: 1.0000e+00 5.3579e+01 1.0616e+02 1.5874e+02 2.1132e+02 Columns 6 through 10: 2.6389e+02 3.1647e+02 3.6905e+02 4.2163e+02 4.7421e+02 Columns 11 through 15: 5.2679e+02 5.7937e+02 6.3195e+02 6.8453e+02 7.3711e+02 Columns 16 through 20: 7.8968e+02 8.4226e+02 8.9484e+02 9.4742e+02 1.0000e+03 octave:33> 1:mrif/m:1000 ans = Columns 1 through 11: 1 51 101 151 201 251 301 351 401 451 501 Columns 12 through 20: 551 601 651 701 751 801 851 901 951 octave:34> catenaria2 octave:35> catenaria2 octave:36> catenaria2 octave:37> catenaria1 octave:38> catenaria1 octave:39> x=[0,1,4,5,7,9,10] x = 0 1 4 5 7 9 10 octave:40> diff(x) ans = 1 3 1 2 2 1 octave:41> h=diff(x) h = 1 3 1 2 2 1 octave:42> [1./h(1);1./h(1:m-2)+1./h(2:m-1);1./h(m-1)] error: A(I): index out of bounds; value 999 out of bound 6 octave:42> m=7 m = 7 octave:43> [1./h(1);1./h(1:m-2)+1./h(2:m-1);1./h(m-1)] error: vertical dimensions mismatch (1x1 vs 1x5) octave:43> h h = 1 3 1 2 2 1 octave:44> x=[0,1,4,5,7,9,10]' x = 0 1 4 5 7 9 10 octave:45> h=diff(x) h = 1 3 1 2 2 1 octave:46> [1./h(1);1./h(1:m-2)+1./h(2:m-1);1./h(m-1)] ans = 1.0000 1.3333 1.3333 1.5000 1.0000 1.5000 1.0000 octave:47> [[-1./h;0],[1./h(1);1./h(1:m-2)+1./h(2:m-1);1./h(m-1)]] ans = -1.00000 1.00000 -0.33333 1.33333 -1.00000 1.33333 -0.50000 1.50000 -0.50000 1.00000 -1.00000 1.50000 0.00000 1.00000 octave:48> x x = 0 1 4 5 7 9 10 octave:49> [[-1./h;0],[1./h(1);1./h(1:m-2)+1./h(2:m-1);1./h(m-1)],[0;-1./h]] ans = -1.00000 1.00000 0.00000 -0.33333 1.33333 -1.00000 -1.00000 1.33333 -0.33333 -0.50000 1.50000 -1.00000 -0.50000 1.00000 -0.50000 -1.00000 1.50000 -0.50000 0.00000 1.00000 -1.00000 octave:50> A=spdiags([[-1./h;0],[1./h(1);1./h(1:m-2)+1./h(2:m-1);1./h(m-1)],[0;-1./h]],[-1,0,1],m,m) A = Compressed Column Sparse (rows = 7, cols = 7, nnz = 19 [39%]) (1, 1) -> 1 (2, 1) -> -1 (1, 2) -> -1 (2, 2) -> 1.3333 (3, 2) -> -0.33333 (2, 3) -> -0.33333 (3, 3) -> 1.3333 (4, 3) -> -1 (3, 4) -> -1 (4, 4) -> 1.5000 (5, 4) -> -0.50000 (4, 5) -> -0.50000 (5, 5) -> 1 (6, 5) -> -0.50000 (5, 6) -> -0.50000 (6, 6) -> 1.5000 (7, 6) -> -1 (6, 7) -> -1 (7, 7) -> 1 octave:51> octave:51> octave:51> octave:51> octave:51> octave:51> octave:51> octave:51> octave:51> octave:51> octave:51> octave:51> full(A) ans = 1.00000 -1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -1.00000 1.33333 -0.33333 0.00000 0.00000 0.00000 0.00000 0.00000 -0.33333 1.33333 -1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -1.00000 1.50000 -0.50000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.50000 1.00000 -0.50000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.50000 1.50000 -1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -1.00000 1.00000 octave:52> A*ones(m,1) ans = 0.0000e+00 -5.5511e-17 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 octave:53> A*[1:m]' ans = -1.00000 0.66667 -0.66667 0.50000 0.00000 -0.50000 1.00000 octave:54> A*x ans = -1 0 0 0 0 0 1 octave:55> diary off