mm-440 === Subject: : Re: minimisation of function with plenty of local minimaYou might also look at the global optimization web site:http://www.mat.univie.ac.at/~neum/glopt.html> folks,>I'm not a expert in numerical analysis and would like to find an efficient>way to minimise a function of a few tens of variables that happens to have>a large number of local minima.>In the Numerical Recipes, I could read sometng about Simulated annealing,>but unfortunately did not grasp a single word of it --- except that it>includes pseudo thermal noise... It also appears that such functions take>a significant number of parameters that I do not underd and that might be>indeed of importance (if not they weren't there I guess).>Could one tell me whether ts type of minimisation is efficient and where>I can read underdable explanation about all ts?=== === Subject: : Re: Quadratical convergence of iterative method> Does anyone know how to check/prove if the iterative method in n-dim space is > quadratically convergent? I've tried using the definition and it doesn't work.> Any idea will be much appreciated.Wch iterative method in n-dimensional space? That phrase describes morethan one algorithm (for more than one problem).=== === Subject: : Re: Quadratical convergence of iterative methodThe method I have looks like ts:x_{k+1} = x_k + f(x_k),where f: R^n --> R^n and is very complicated. Ceculy> Wch iterative method in n-dimensional space? That phrase describes more> than one algorithm (for more than one problem).> ---=== === Subject: : Re: Quadratical convergence of iterative method >The method I have looks like ts:x_{k+1} = x_k + f(x_k),where f: R^n --> R^n and is very complicated. >a _sufficient_ conditoion for quadratic convergence in ts case is: f is two times differentiable in some neighborhood of a zero x*(wch will be a fixed point of the iteration) and the Jacobian (first derivative)is zero there.apply Taylors theorem. === === Subject: : Iterative solution of Ax=b, what's wrong?I have a system of linear equations, Ax=b. I know that A isnon-singular, and I have also solved the system with gaussianelimination, just to make sure.Actually, what I want is to solve it iteratively. I rewrite the systemas follows:x_{k+1} = -inv(diag(A)) * (A - diag(A)) * x_k + inv(diag(A)) * bThe algorithm converges, but not to the solution. What's wrong?Round-off? Could it really be that lim x_k is NOT the solution ofAx=b?In my example, A is a 6x6 matrix and its elements range from 10^-5 to10^2.I do not tnk that I've made any implementation mistakes. When Isubstitute the solution obtained with gaussian elimination in theiterative formula above, I get what I expect. for any nts.=== === Subject: : Re: Iterative solution of Ax=b, what's wrong?I have a system of linear equations, Ax=b. I know that A is> non-singular, and I have also solved the system with gaussian> elimination, just to make sure.Actually, what I want is to solve it iteratively. I rewrite the system> as follows:x_{k+1} = -inv(diag(A)) * (A - diag(A)) * x_k + inv(diag(A)) * bThe algorithm converges, but not to the solution. What's wrong?> Round-off? Could it really be that lim x_k is NOT the solution of> Ax=b?In my example, A is a 6x6 matrix and its elements range from 10^-5 to> 10^2.I do not tnk that I've made any implementation mistakes. When I> substitute the solution obtained with gaussian elimination in the> iterative formula above, I get what I expect. for any nts.> Sincex_{k+1} - x_k = inv(diag(A)) * (b - A * x_k)if it converges then it does so to the correct answer!The most likely explanation is a trivial but difficult to spot error.You'd probably be best posting the code. I'd volunteer to have a lookbut then it would turn out to be implemented in some language of wchI've never heard.=== === Subject: : Re: Iterative solution of Ax=b, what's wrong?> Since> x_{k+1} - x_k = inv(diag(A)) * (b - A * x_k)> if it converges then it does so to the correct answer!> The most likely explanation is a trivial but difficult to spot error.> You'd probably be best posting the code. I'd volunteer to have a look> but then it would turn out to be implemented in some language of wch> I've never heard.> for your answer and offer. There was indeed an implementationerror wch I've fixed now.Still, I have a more general question:How does one solve Ax=b iteratively in the general case ('A' non-singular,non-symetric and maybe not positive definite)? Does a relaxation parameterw always exist such that the SOR algorithm converges? If not, wchiterative algorithms work in the general case? in advance,=== === Subject: : Re: Iterative solution of Ax=b, what's wrong?Since x_{k+1} - x_k = inv(diag(A)) * (b - A * x_k) if it converges then it does so to the correct answer! The most likely explanation is a trivial but difficult to spot error.> You'd probably be best posting the code. I'd volunteer to have a look> but then it would turn out to be implemented in some language of wch> I've never heard. for your answer and offer. There was indeed an implementation> error wch I've fixed now.Still, I have a more general question:How does one solve Ax=b iteratively in the general case ('A' non-singular,> non-symetric and maybe not positive definite)? Does a relaxation parameter> w always exist such that the SOR algorithm converges? If not, wch> iterative algorithms work in the general case? in advance,> The short answer is that a relaxation parameter w does not alwaysexist such that SOR versions of Jacobi or Gauss-Seidel algorithmsconverge.However all is not doom and gloom. Suppose we have an initial estimateof the solution x1 and an estimate x2 generated by our iterativemethod (whatever it may be - Jacobi, Gauss-Seidel...).Now choose y3 = l1x1 + l2x2 where l1 and l2 are chosen to minimise (b l1Ax1 l2Ax2)T(b l1Ax1 l2Ax2) (T is supposed to be thetranspose operator).i.e. l1 & l2 satisfy the equationsx1ATAx1. l1 + x1ATAx2. l2 = bTAx1x1ATAx2. l1 + x2ATAx2. l2 = bTAx2equations above tox1ATAx1. l1 + x1ATAx2. l2 + x1ATAx3. l3 = bTAx1x1ATAx2. l1 + x2ATAx2. l2 + x2ATAx3. l3 = bTAx2x1ATAx3. l1 + x2ATAx3. l2 + x3ATAx3. l3 = bTAx3(assuming the x3 generated by our iterative method is linearlyindependent of x1 and x2) ...Now the sequence x1,x2,x3... will converge to the correct answerunless the sequence of x's becomes linearly dependent. If A isnon-singular then the iterative method used would have to be prettypoor, for example if it always selected the same value. The reason itwill converge is simple enough - eventually you will have solved anequivalent set of equations.Ts is not really a practical method as it ds. In practice xnwill tend to be linearly dependent on previous values and theequations will become singular. Also there's more and more workinvolved in generating the next set of equations. However, methodsusing the last 2 or 3 estimates don't require too much effort and cansignificantly improve the rate of convergence, particularly in caseswhere convergence is poor.I have used a simple version of it involving repeated solutions of 2x2equations to force the Gauss-Seidel method to converge when itotherwise would not. Ts is effectively the same as SOR but withdifferent relaxation parameters chosen for each iteration.=== === Subject: : Re: Iterative solution of Ax=b, what's wrong?> How does one solve Ax=b iteratively in the general case ('A' non-singular,> non-symetric and maybe not positive definite)? Does a relaxation parameter> w always exist such that the SOR algorithm converges? If not, wch> iterative algorithms work in the general case?You've already received a pointer to the templates book.The almost-trivial answer is that A^t Ax = A^t b has a symmetricpositive definite coefficient matrix, so Conjugate Gradients -- and ahost of related methods -- should work. In exact arithmetic. And dontask at what convergence rate. The general answer is that ts is a hard question. Most iterativemethods depend in some way on eigenvalues of the matrix. If your matrixis really badly conditioned, you can try first to temper it with apreconditioner, but that requires that the matrix is amenable toreasonable approximation to begin with. In some cases it's just veryhard to apply iterative methods. For ince, CG-like methods are basedon the fact that with orthogonality you get projection and thereforeminimisation (imagine me waving my hands very hard wle I type ts),but that means that your matrix has to be at least a bit like elliptic(and no, I'm not going to be any more precise here). Now consider aconvection-diffusion equation, wch has itty-bitty ellipticity andbunches of skew-symmetry. Those are hard. Even worse, take large systemsof ODEs, with a predictor-corrector or so. Those matrices are notng atall like elliptic, so iterative methods work abysmally.I'm sure ts was of no help to you. === === Subject: : Re: Iterative solution of Ax=b, what's wrong?> How does one solve Ax=b iteratively in the general case ('A' non-singular,> non-symetric and maybe not positive definite)? Does a relaxation parameter> w always exist such that the SOR algorithm converges? If not, wch> iterative algorithms work in the general case?If iterative methods interest you, a nice overview can be found at Possible candidate solvers are GMRES or QMR, for your kind of matrix.-- === === Subject: : Re: Iterative solution of Ax=b, what's wrong?I have a system of linear equations, Ax=b. I know that A is >non-singular, and I have also solved the system with gaussian >elimination, just to make sure.Actually, what I want is to solve it iteratively. I rewrite the system >as follows:x_{k+1} = -inv(diag(A)) * (A - diag(A)) * x_k + inv(diag(A)) * b jacobi iteration The algorithm converges, but not to the solution. What's wrong? >Round-off? Could it really be that lim x_k is NOT the solution of >Ax=b?In my example, A is a 6x6 matrix and its elements range from 10^-5 to >10^2.I do not tnk that I've made any implementation mistakes. When I >substitute the solution obtained with gaussian elimination in the >iterative formula above, I get what I expect. >assume that it converges to x*. then x*= -inv(diag(A)) * (A - diag(A)) * x* + inv(diag(A)) * bmultiply by diag(A), cancel diag(A)*x* on both sieds then 0=-Ax*+bhence x* is a solution. hence either the system is solvable but not uniquelysolvable hence the matrix singualr and since you state it is not then x* must be the solution. but what is convergence numerically for you.ts iteration might be so slow that you tnk to be at the limit being farfrom it. hence analyze the iegenvalues of inv(diag(A)) * (A - diag(A))wch determine the speed of convergence .=== === Subject: : Solve Ax=b if A and b are relying on x I encounter an optimization problem wch I converted to a linear system form like Ax=b.But here Both A and b are relying on x. Youcan consider they as function of (x), i.e. A(x)x = b(x). However, there is no close formfor either A(x) and b(x). Anyway, we can assumethat A(x) and b(x) are C^2. Can we still apply any iterative methodto solve such problem? For example, startingfrom x0, after one iteration according toA0 and b0, get x1, and then update A and baccordingly to get next x2 and so on? OrI have to use nonlinear optimization method?=== === Subject: : Re: Solve Ax=b if A and b are relying on x > I encounter an optimization problem wch >I converted to a linear system form like Ax=b. >But here Both A and b are relying on x. You >can consider they as function of (x), i.e. >A(x)x = b(x). However, there is no close form >for either A(x) and b(x). Anyway, we can assume >that A(x) and b(x) are C^2. > Can we still apply any iterative method >to solve such problem? For example, starting >from x0, after one iteration according to >A0 and b0, get x1, and then update A and b >accordingly to get next x2 and so on? Or >I have to use nonlinear optimization method?your ideas lead to the iteration x(k+1) = inverse(A(x(k))*b(x(k)) def=F(x(k)) and the convergence of ts iteration can be investigated by considering the Jacobian matrix of the right hand side with respect to x. If at the solutionall its eigenvalues have absolute value less than one, then the iterationconverges locally (in a small neighborhood of the solution) at least.ts is Ostrowkis theorem. The Jacobian has as i-th column (partial/partial x_i ) F(x) = -inverse(A(x))*{(partial/partial x_i)A(x)}*inverse(A(x))*b(x)+ inverse(A(x))*( (partial/partial x_i) b(x) )in most cases you will see that Ostrowkis theorem is not applicable. you could still try to solve the nonlinear system G(x)=A(x)*x-b(x) = 0 by Newton's method or some of its variants, but again ts is not a good idea if you have no good initial guess. Since you write that the system comes from an optimization problem I guess these are the necessary optimality conditions.using these directly for solving the problem is one of the worst ideas one canhave. there is plenty of good software for solving linear and nonlinear optimizationproblems, see e.g.http://plato.la.asu.edu/guide.html === === Subject: : recurrence relation evaluation trouble,I have a three-term recurrence relation. The solution I want is the set oforthogonal polynomials originating from ts recurrence relation(normalisation : p_0 = 1). I know the weight function and the interval oforthogonality. Orthogonality has been proved. First issue: the interval is [-1,+1] and two discrete points (pm nu_0)where the weight function is a Dirac impulse. The position of these twopoints is mainly governed by a parameter c in [0,1]. The closer c to 1,the larger |nu_0|.Second issue: I have no explicit form for the polynomials, only therecurrence relation for evaluation.Trd issue: I need to evaluate these polynomials in these two discretepoints but these discrete points become very good approximations for thelargest root of the polynomials (because the dice between the largestroot of an orthogonal polynomial and the right edge point of theorthogonality interval decreases with n). Fourth issue: the recurrence relation contains a term with (1-c). For cclose to one (I'm talking about c = 0.99999; 1-c has only 12 correct digitsin double precision). Combined with the fact that I'm evaluating very closeto a root, ts gives cancellation. I implemented the recurrence relation (forward)and already for the 3rddegree polynomial I have only 5 significant digits left, the 4th degreepolynomial is completely bogus. When I use a quadruple double precisionpackage, I get the 4th degree with 6 digits, but the 5th degree is againbogus. I tried the backward approach (Miller) but the problem is that thecoefficients in the recurrence don't behave like good asymptotics and tsgives very poor results.Can somebody advise me on a more stable evaluation technique?=== === Subject: : Re: recurrence relation evaluation trouble >, >I have a three-term recurrence relation. The solution I want is the set of >orthogonal polynomials originating from ts recurrence relation >(normalisation : p_0 = 1). I know the weight function and the interval of >orthogonality. Orthogonality has been proved. First issue: the interval is [-1,+1] and two discrete points (pm nu_0) >where the weight function is a Dirac impulse. The position of these two >points is mainly governed by a parameter c in [0,1]. The closer c to 1, >the larger |nu_0|. Second issue: I have no explicit form for the polynomials, only the >recurrence relation for evaluation.Trd issue: I need to evaluate these polynomials in these two discrete >points but these discrete points become very good approximations for the >largest root of the polynomials (because the dice between the largest >root of an orthogonal polynomial and the right edge point of the >orthogonality interval decreases with n). Fourth issue: the recurrence relation contains a term with (1-c). For c > close to one (I'm talking about c = 0.99999; 1-c has only 12 correct digits >in double precision). Combined with the fact that I'm evaluating very close >to a root, ts gives cancellation. I implemented the recurrence relation (forward)and already for the 3rd >degree polynomial I have only 5 significant digits left, the 4th degree >polynomial is completely bogus. When I use a quadruple double precision >package, I get the 4th degree with 6 digits, but the 5th degree is again >bogus. I tried the backward approach (Miller) but the problem is that the >coefficients in the recurrence don't behave like good asymptotics and ts >gives very poor results.Can somebody advise me on a more stable evaluation technique? ??????? your scalar product is (f,g)=f(-ny)*g(-ny)+f(ny)*g(ny) ts would give only two independent polynomials? the set of zeroes of the weight function must have measure zero on [-1,1] in your case. but if the weight function behaves only approximately like a Dirac impulse then you have indeed a problem with the numerics. Walter Gautsc did a lot of work on computing orthogonal polynomials in a stable manner (clearly _not_ using the recursion) and s software Orthpol in http://www,netlib.org/toms/726 might help you. but the case is indeed a hard one then and might nevertheless require an extreme precision in the computation. === === Subject: : Re: recurrence relation evaluation trouble >, I have a three-term recurrence relation. The solution I want is the set oforthogonal polynomials originating from ts recurrence relation(normalisation : p_0 = 1). I know the weight function and the interval oforthogonality. Orthogonality has been proved.First issue: the interval is [-1,+1] and two discrete points (pm nu_0)where the weight function is a Dirac impulse. The position of these twopoints is mainly governed by a parameter c in [0,1]. The closer c to 1,the larger |nu_0|. Second issue: I have no explicit form for the polynomials, only therecurrence relation for evaluation.Trd issue: I need to evaluate these polynomials in these two discretepoints but these discrete points become very good approximations for thelargest root of the polynomials (because the dice between the largestroot of an orthogonal polynomial and the right edge point of theorthogonality interval decreases with n).Fourth issue: the recurrence relation contains a term with (1-c). For c close to one (I'm talking about c = 0.99999; 1-c has only 12 correct digitsin double precision). Combined with the fact that I'm evaluating very closeto a root, ts gives cancellation.I implemented the recurrence relation (forward)and already for the 3rddegree polynomial I have only 5 significant digits left, the 4th degreepolynomial is completely bogus. When I use a quadruple double precisionpackage, I get the 4th degree with 6 digits, but the 5th degree is againbogus. I tried the backward approach (Miller) but the problem is that thecoefficients in the recurrence don't behave like good asymptotics and tsgives very poor results.Can somebody advise me on a more stable evaluation technique? ??????? your scalar product is> (f,g)=f(-ny)*g(-ny)+f(ny)*g(ny)> ts would give only two independent polynomials?> the set of zeroes of the weight function must have measure zero on [-1,1] in your> case. but if the weight function behaves only approximately like a Dirac impulse> then you have indeed a problem with the numerics. Walter Gautsc did a lot> of work on computing orthogonal polynomials in a stable manner> (clearly _not_ using the recursion) and s software Orthpol> in> http://www,netlib.org/toms/726> might help you. but the case is indeed a hard one then and might nevertheless> require an extreme precision in the computation. > If you need gh precision, check out the recent column in Computingin Science and Engineering (CiSE). Author was David Smith.-- ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/=== === Subject: : Re: recurrence relation evaluation trouble> ??????? your scalar product is> (f,g)=f(-ny)*g(-ny)+f(ny)*g(ny)> ts would give only two independent polynomials?no. the scaler product is int(w*f*g,x=-1..1) + M1*f(nu0)*g(nu0) +M2*f(-nu0)*g(-nu0). You can see ts as an integral on the interval[-nu0,nu0] with a weight function that consists of the function w wch isonly non-zero in the sub-interval [-1,1] and two Dirac impulses at pm nu0.> the set of zeroes of the weight function must have measure zero on [-1,1]> in your case. but if the weight function behaves only approximately like> a Dirac impulse then you have indeed a problem with the numerics. Walter> Gautsc did a lot of work on computing orthogonal polynomials in a> stable manner (clearly _not_ using the recursion) and s software > http://www,netlib.org/toms/726=== === Subject: : Pseudoinvers Matrix Computation Implemented in C sourceI'm writing a statistical appliction wch requires the (possiblyfast) computation of the pseudoinverse matrix of a 6x500 matrix ofreal numbers:sometimes the results are incorrect (I'm comparing the results to myoriginal algorithm, written in Matlab), and i found ts is due toclose-to-singularity problems, wch are avoided in Matlab enviromentby their implementation of the pseudoinverse computation:computing optimization or scaling, required to fix theC-implementation problem; or, even better, a C function computingpseudiinverse that I can use or study?Filippo Venturi=== === Subject: : Re: Pseudoinvers Matrix Computation Implemented in C sourceI'm writing a statistical appliction wch requires the (possibly >fast) computation of the pseudoinverse matrix of a 6x500 matrix of >real numbers: >sometimes the results are incorrect (I'm comparing the results to my >original algorithm, written in Matlab), and i found ts is due to >close-to-singularity problems, wch are avoided in Matlab enviroment >by their implementation of the pseudoinverse computation: >computing optimization or scaling, required to fix the >C-implementation problem; or, even better, a C function computing >pseudiinverse that I can use or study? >Filippo Venturiclapack has a svd. (c tranlation from lapack)http://www.netlib.org/clapack/double/dgelss.csince your problem is small, ts is all you need.=== === Subject: : Re: Pseudoinvers Matrix Computation Implemented in C source>I'm writing a statistical appliction wch requires the (possibly>fast) computation of the pseudoinverse matrix of a 6x500 matrix of>real numbers:>sometimes the results are incorrect (I'm comparing the results to my>original algorithm, written in Matlab), and i found ts is due to>close-to-singularity problems, wch are avoided in Matlab enviroment>by their implementation of the pseudoinverse computation:>computing optimization or scaling, required to fix the>C-implementation problem; or, even better, a C function computing>pseudiinverse that I can use or study?>Filippo VenturiI believe people usually use sometng like QR factorisation to do tssort of tng.-- === === Subject: : matlab programming questionI can't figure out how to program the following in matlab:N = positive integerH(1,1, ..., 1).matrix = 1 % (H has N arguments)The problem is that you need N for-loops. You should my your own loop-function somehow, but I don't know how?Wilbert=== === Subject: : Re: matlab programming question> I can't figure out how to program the following > in matlab:N = positive integer> H(1,1, ..., 1).matrix = 1 % (H has N arguments)The problem is that you need N for-loops. You > should my your own loop-function somehow, but I > don't know how?How exactly did you define H in matlab? If I underdts right, you want H to be sometng like an arrayof matrices? Or is H an nxm matrix and you want to setall elements with n=1 to 0?maybe somebody in comp.soft-sys.matlab can help you...Nils=== === Subject: : Re: matlab programming question> I can't figure out how to program the following> in matlab: N = positive integer> H(1,1, ..., 1).matrix = 1 % (H has N arguments) The problem is that you need N for-loops. You> should my your own loop-function somehow, but I> don't know how?How exactly did you define H in matlab? If I underd> ts right, you want H to be sometng like an array> of matrices? Or is H an nxm matrix and you want to set> all elements with n=1 to 0?maybe somebody in comp.soft-sys.matlab can help you...yet, maybe it should be a structure.What I want is the following:[y] = function(N) % N integerH = H(a_j1,a_j2, ..., a_jN) % jk runs from 1 to Ny(a_j1,a_j2, ..., a_jN) = sqrt(a_j1^2 + a_j2^2 + ... + a_jN^2)% thus y has N^N values, and must be a NxNx...xN (N times) matrix.Wilbert=== === Subject: : Re: matlab programming questionThe problem is that 3 is not arbitrary. What I want is to make the> following function (or one wch produces the same output):[h] = function(N)> for j1=1:3, ..., for jN=1:3,> h(j1, ..., jN) = (j1 + ... + jN)*sqrt(j1^2 + ... + jN^2);> end; ... end;How about the following implementation:function h = hfunc(N,M)% N is number of dimensions.% M is max index along each dimension (3 in your example).index = 1:M; % Indices take on these valuesxstr = '[x1'; % Initializefor i = 2:N % Loop over dimensions xstr = sprintf('%s,x%d', xstr, i); % Build up LHS stringendxstr = [xstr ']']; % Terminate the string properly.eval([xstr ' = ndgrid(index);']); % Create N matrices.sum1 = zeros(size(x1)); % Initializationsum2 = zeros(size(x1)); % Initializationfor i = 1:N eval(sprintf('sum1 = sum1 + x%d;', i)); eval(sprintf('sum2 = sum2 + x%d.^2;', i));endh = sqrt(sum2) .* sum1;returnTs produces the following output:> h = hfunc(2,2)h = 2.8284 6.7082 6.7082 11.3137> h = hfunc(2,3)h = 2.8284 6.7082 12.6491 6.7082 11.3137 18.0278 12.6491 18.0278 25.4558> h = hfunc(3,3)h(:,:,1) = 5.1962 9.7980 16.5831 9.7980 15.0000 22.4499 16.5831 22.4499 30.5123h(:,:,2) = 9.7980 15.0000 22.4499 15.0000 20.7846 28.8617 22.4499 28.8617 37.5233h(:,:,3) = 16.5831 22.4499 30.5123 22.4499 28.8617 37.5233 30.5123 37.5233 46.7654Hope ts helps,-- Simon=== === Subject: : Re: matlab programming question> What I want is the following:> [y] = function(N) % N integer> H = H(a_j1,a_j2, ..., a_jN)> % jk runs from 1 to N> y(a_j1,a_j2, ..., a_jN) = sqrt(a_j1^2 + a_j2^2 + ... + a_jN^2)Wilbert, I can't say that I completely underd what you're trying to do,but it seems like you are asking how to index into H with N indices, when Nis not known until run time. You can use a cell array for that. Here's anexample where I create a 3x3x3 H array, but then I index into it withoutexplicitly writing an expression having 3 indices:> for i=1:3, for j=1:3, for k=1:3, h(i,j,k) = 100*i + 10*j + k; end; end;end> ind = {2,1,3};> h(ind{:})ans = 213If that's not the gist of your question, please let me know.-- Tom=== === Subject: : Re: matlab programming questionWhat I want is the following: [y] = function(N) % N integer H = H(a_j1,a_j2, ..., a_jN)> % jk runs from 1 to N y(a_j1,a_j2, ..., a_jN) = sqrt(a_j1^2 + a_j2^2 + ... + a_jN^2)Wilbert, I can't say that I completely underd what you're trying to do,> but it seems like you are asking how to index into H with N indices, when N> is not known until run time. You can use a cell array for that. Here's an> example where I create a 3x3x3 H array, but then I index into it without> explicitly writing an expression having 3 indices: for i=1:3, for j=1:3, for k=1:3, h(i,j,k) = 100*i + 10*j + k; end; end;> end> ind = {2,1,3};> h(ind{:})> ans => 213If that's not the gist of your question, please let me know.The problem is that 3 is not arbitrary. What I want is to make the following function (or one wch produces the same output):[h] = function(N)for j1=1:3, ..., for jN=1:3, h(j1, ..., jN) = (j1 + ... + jN)*sqrt(j1^2 + ... + jN^2); end; ... end;I hope it's more clear now. I don't see how to use your nt above forsolving ts.Wilbert=== === Subject: : Re: matlab programming question> The problem is that 3 is not arbitrary. What I want is to make the> following function (or one wch produces the same output):> [h] = function(N)> for j1=1:3, ..., for jN=1:3,> h(j1, ..., jN) = (j1 + ... + jN)*sqrt(j1^2 + ... + jN^2);> end; ... end;Okay, now I'm intrigued. Here's a function that computes the function youwant, and assigns it to output arrays two different ways. Let me know ifthat doesn't get at the heart of your question.-- Tomfunction [h1,h2]=doit(N)% Do Wilbert's tngbase = 3.^(N-1:-1:0); % base for converting to indicesh1 = zeros(repmat(3,1,N)); % initialize output to 0h2 = h1;for k=0:3^N-1 ind = 1+mod(floor(k./base),3); % get indices result = sum(ind) * sqrt(sum(ind.^2)); % compute Wilbert's function h2(k+1) = result; % put into array w/ single index c = num2cell(ind); % or get cell array of indices h1(c{:}) = result; % assign with multiple indicesend=== === Subject: : Re: matlab programming questionquas litteras Wilbert Dijkhof in continuo sci.math.num-analysis:66467 scripsit eis responsurus sum> N = positive integer> H(1,1, ..., 1).matrix = 1 % (H has N arguments)The problem is that you need N for-loops. I do not know matlab, but ts a a more general problem.If you software allow 1-D-indexing of multidimensional quantities it'sstraighforward -- but you have to determine whether you deal with arow-major data storing or a column-major one. Is a(i,j) the sameas a(i+a*j) or as a(a*i+j)?If not, you can do the following with only two loops index = [0, 0, ...]; for (int i = 1; i < k^N; ++i) { index(1) += 1; for (int j = 1; j < N; ++j) { if (index(j) == k+1) { index(j+1) += 1; index(j) = 0; } } H (index(1), ..., index(n)).matrix = whatever; }[index can be seen as a base k number to wch you add 1 at each iteration.]Btw, what's the link between your question, sauerkraut, and numericalanalysis?-- R.8egis=== === Subject: : Suggested Solver for Diffusion Type ProblemHelloI need to solve a large set of equations wch are obtained by applying afinite difference formulation to a rectangular grid - actually a greyscaleimage. The equations are based on minimising the magnitude of a Laplacianoperator and are sparse and diagonally dominant.So far I've been using a simple Conjugate Gradient solver without anypre-conditioning. Ts works but is very slow. To speed it up I've taken amultigrid approach - reducing the resolution to 8x8 and solving atprogressively gher resolutions. Ts works much better but is still a bitslow at the gher resolutions.Since ts is a fairly dard problem I was wondering if there is a wellestablished solver wch is optimised for ts type of problem?Pete=== === Subject: : Signal synthesis from Wigner-Ville distributionI have some problems with ts question. Can anybody help me?I'm trying to recover signal from Wigner-Ville distribution aftertime-frequency filtering. I base on Spatial and Time-FrequencySignature Estimation of Nonstationary Sources paper [Moeness] (foundin internet).The procedures from ts paper work as follows:1. Calculate Wigner-Ville distribution; WVD(t,f) = sum[ x(t+k/2)*x'(t-k/2)*exp(-jkf) ]2. Filter WVD as you need WVD'(t,f) = WVD(t,f)*G(t,f) where G(t,f) - binary mask (0,1) 3. Take the inverse Fast Fourier Transform of WVD' P(t,n) = IFFT( WVD' ) 4. Construct the matrix Q with Q(i,j) = P( (i+j)/2 , i-j )5. Apply eigen-decomposition to the matrix Q and obtain the maximumeigenvalue and associated eigenvector. Now recovered signal should be based on ts eigenvector. And now my computation: i.e. signal: sig=[1 2 3 4 5 6 7 1] (WVD should obtained from signal transformed by lbert Transform,but I want to make it easier) WVD = Columns 1 through 7 1.0000 10.0000 35.0000 84.0000 119.0000 114.0000 61.0000 1.0000 8.2426 20.3137 27.3137 56.1127 85.4975 57.4853 1.0000 4.0000 -1.0000 -8.0000 -17.0000 28.0000 49.0000 1.0000 -0.2426 -2.3137 4.6863 -6.1127 -13.4975 40.5147 1.0000 -2.0000 3.0000 -4.0000 15.0000 -26.0000 37.0000 1.0000 -0.2426 -2.3137 4.6863 -6.1127 -13.4975 40.5147 1.0000 4.0000 -1.0000 -8.0000 -17.0000 28.0000 49.0000 1.0000 8.2426 20.3137 27.3137 56.1127 85.4975 57.4853 Column 8 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00002. I'm going to do ts without mask to make it easier3. My inwerse transform P= Columns 1 through 7 1.0000 4.0000 9.0000 16.0000 25.0000 36.0000 49.0000 0 3.0000 8.0000 15.0000 24.0000 35.0000 6.0000 0 0 5.0000 12.0000 21.0000 4.0000 0 0 0 0.0000 7.0000 2.0000 0 0.0000 0 0 0 0 0 0 0 0 0 0.0000 7.0000 2.0000 0 0.0000 0 0 5.0000 12.0000 21.0000 4.0000 0 0 3.0000 8.0000 15.0000 24.0000 35.0000 6.0000 Column 8 1.0000 0 0 0 0 0 0 04. And in the end Q-matrix Q= Columns 1 through 7 2.0000 0 3.0000 0 5.0000 0 7.0000 0 8.0000 0 8.0000 0 12.0000 0 3.0000 0 18.0000 0 15.0000 0 21.0000 0 8.0000 0 32.0000 0 24.0000 0 5.0000 0 15.0000 0 50.0000 0 35.0000 0 12.0000 0 24.0000 0 72.0000 0 7.0000 0 21.0000 0 35.0000 0 98.0000 0 2.0000 0 4.0000 0 6.0000 0 Column 8 0 2.0000 0 4.0000 0 6.0000 0 2.0000 5. After that eigenvectors/eigenvalues are copmputed by Matlab, soit's no problem The results obtained as above are not satisfied for me because whenI try to use filtering the recovered signal is not correct in timedomain.I tnk WVD and P matrix are computed correctly but I'm not sure aboutQ-matrix.Has anybody use ts method and can help to solve my problem?Maybe you know another method that work in signal synthesis fromWigner-Ville distribution?I'll be grateful for any help. === === Subject: : The Analytic Continuation of the HyperPower FunctionThe LambertW can be used to solve a couple of interesting problems,including some relating to the fixed points of the hyperpower sequencex^x^x^...Turns out there is another bonus in there, that's again using theLambertW, ts time to define the analytic continuation of the realhyperpower function F(x) = x^x^x^...investigation, but the idea that it could be used as the analyticcontinuation of F(x) hadn't occured to me, until two days ago.There is also a nice investigation of the two components of the Hopfbifurcation using Maple, and some code to sketch the branches in(0,e^(-e)),in my math section:Enjoy.-- Ioannishttp://users.forthnet.gr/ath/jgal/ === === Subject: : S-Lang modules for GNU Scientific Library (GSL) released I am pleased to announce the availability of a set of S-Langmodules for the GNU Scientific Library. A brief introduction to themodules is attached below. For additional information about theiruse and availability, them, please visit . The GNU Scientific Library (GSL ) is a vast collection of robust and well documented numerical functions. It includes support for many special functions, random numbers, interpolation and integration routines, and much more. For more information about GSL, visit . Many of the routines in the GSL may be made available to the S-lang interpreter via the GSL modules described by ts document, whose most recent version may be found at . At the moment, four GSL modules are available: o gslsf: The GSL special function module. Currently, ts module provides an interface to nearly 200 GSL special functions. o gslconst: The GSL conts module. Ts module defines many conts such as CONST_MKS_SPEED_OF_LIGHT, CONST_CGS_BOLTZMANN, etc. o gslinterp: The GSL interpolation module, wch includes routines for linear interpolation, cubic splines, etc. o gslcore: Ts is a module that must be loaded before any of the above modules can be loaded. Its main purpose is to provide support functions for the other GSL modules. There are many functions that are not yet wrapped. For example, none of GSL's random number routines have been wrapped. Future releases of the GSL module will include more functionality. Nevertheless, what has been implemented should prove useful. -- _______________________________________________To unsubscribe, visit http://jedsoft.org/slang/mailinglists.html=== === Subject: : Testing numerical softwareDistribution: worldDoes somebody recommend journal artical or books on the designof test cases for numerical software?Christian -- Dipl.-Ing. Christian KeckPhysikalisch-Technische BundesanstaltFachlaboratorium 5.32 (Koordinatenmessgeraete)Bundesallee 10038116 BraunschweigTelefon (+49) (0531) 592-5323=== === Subject: : Re: Testing numerical software> Does somebody recommend journal artical or books on the design> of test cases for numerical software?You may find it helpful to search the on-line bookstores on verificationand validation.=== === Subject: : Re: Testing numerical softwareDistribution: world>Does somebody recommend journal artical or books on the design>of test cases for numerical software?>Christian>-- a search in zentralblatt fuer mathmeatik with testing and softwaregave 79 nts, notng about numerical software. the reason is thatspecial software needs special testing techniques. there is a lot a discussion about benching e.g. in optimization (e.g. J.J . Mreo`where are the Handbooks of Westlake and Gregory & Karney for numericlalinear algebra testcases, there are the collections of testcases assembled for optimization e.g. by Floudas et al (a book), (see http://plato.la.au.edu/topics/benchm.html)there are collections of differential equations problems for testing ode software, but notng general. if you look in ZBL nder numerical and software you get 82 ts,alomost anytng special.sometng genral: 7. Zbl 0866.65001 Ueberhuber, Christoph W.Numerical computation. Methods, software, and analysis. In 2 vol. (English)Berlin: Springer; 3-540-62057-5 3. Zbl 0932.00014 Redivo-Zaglia, M.(ed.)Numerical software. (English)Numer. Algorithms 20, No.2-3, 107-268 (1999). MSC 2000: *00B15 65-0619. Zbl 0863.00021 Daehlen, Morten (ed.); Tveito, Aslak (ed.)Numerical methods and software tools in industrial mathematics. (English)Boston, MA: Birkh.8auser. 400 p. DM 148,00; .9aS 1.081,00; sFr 128,00 (1997). MSC 2000: *00B15 65-0620. Zbl pre01069521 Boyle, James M.; Harmer, Terence J.; Winter, Victor L.Arge, Erlend (ed.) et al., Modern software tools for scientific computing. International workshop, Oslo, Norway, September 16-18, 1996. Boston: Birkh.8auser. 353-372 (1997). MSC 2000: *68-9935. Zbl 0791.68031 Lucena, Carlos J.P.; Qian, YueminFormal specification methods and numerical software: A case study using $Z$. (English)Util. Math. 44, 85-114 (1993). MSC 2000: *68N99 68Q60 65F506. Zbl 0782.65004 Rice, R.Numerical methods, software, and analysis.2. ed. (English)Boston, MA: Academic Press, Inc.. xiv, 720 p. (1993). MSC 2000: *65-01 65-04, Reviewer: V.Mehrmann42. Zbl 0825.68261 Jacobs, D.A.H.; ham, G.Experiences with some software engineering practices in numerical software. (English)Reliable numerical computation, Proc. Conf. in Honour of J. H. Wilkinson, Teddington/UK, 277-296 (1990). MSC 2000: *68N01 65Y9944. Zbl 0718.65013 Dongarra, Jack; Hammarling, SvenEvolution of numerical software for dense linear algebra. (English)Reliable numerical computation, Proc. Conf. in Honour of J. H. Wilkinson, Teddington/UK, 297-327 (1990). MSC 2000: *65F05 65-03 65Y05, Reviewer: H.R.Schwarz=== === Subject: : recursive loops and memoryhello group I had a question about loop constructs and random numbergenerators(I'm using mathematica but th eidea shoudl hold for anynumerics applications)PLease consider following construct.For[it=1, it<= 10, NotebookCreate[]; For[ii=1, ii<=100, Label[sample] pick Random[RealNumbers]; use Random Real numbers to plug into ODE system; Check[ for errors made by ODE solver and plotters, IF errors, go to Label[sample] Clear[Random numbers and interpolation madeby NDSolve] IF no errors[plot and output to another notebook] ] ]Above construct runs for a single outer For Loop iteration, then ithangs.... and then it runs out of memory...is it because of the inner recursive loop of picking random realnumbers?Ts post really belongs in Mathematica forum but the usenet serverseems to be down for sometime now.any comments or suggestions is thoroughly appreciated. sean=== === Subject: : Re: recursive loops and memoryhello groupI had a question about loop constructs and random number> generators(I'm using mathematica but th eidea shoudl hold for any> numerics applications)PLease consider following construct.For[it=1, it<= 10,> NotebookCreate[]; For[ii=1, ii<=100,> Label[sample]> pick Random[RealNumbers];> use Random Real numbers to plug into ODE system;> Check[ for errors made by ODE solver and plotters,> IF errors, go to Label[sample]> Clear[Random numbers and interpolation made> by NDSolve]> IF no errors[plot and output to another notebook]> ] ]Above construct runs for a single outer For Loop iteration, then it> hangs.... and then it runs out of memory...is it because of the inner recursive loop of picking random real> numbers?Ts post really belongs in Mathematica forum but the usenet server> seems to be down for sometime now.any comments or suggestions is thoroughly appreciated.seanRecursion is implemented on a computer using the stack. You canrun out of stack if the recursion gets too deeply nested.I suggest a look at my column, Recurses! that appeared in Computingin Science and Engineering (IEEE/AIP).You will find it athttp://www.phys.virginia.edu/classes/551.jvn.fall01/CiSE_ progs/Cprogs.htm-- ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/