mm-474 === Subject: Measure for oscillating time seriesI'm trying to find a suitable function to measure the jaggedness of aline. The obvious idea seems to be to use the standard deviation, butthis is not giving me the results I require.Essentially: - If the series oscillates by more than 1 point many times themeasure should increase. - If the series makes large smooth oscillations (e.g. goes instraight lines from -5 to 5 and back again) the measure should godown. - If the series makes large but jagged oscillations (e.g.5*sin(x/30)+0.1*sin(x) then the measure should be higher than theabove.Could anyone here help me to find a solution to this?Anon. === Subject: Re: Measure for oscillating time series> I'm trying to find a suitable function to measure the jaggedness of a> line. The obvious idea seems to be to use the standard deviation, but> this is not giving me the results I require.> Essentially:> - If the series oscillates by more than 1 point many times the> measure should increase.> - If the series makes large smooth oscillations (e.g. goes in> straight lines from -5 to 5 and back again) the measure should go> down.> - If the series makes large but jagged oscillations (e.g.> 5*sin(x/30)+0.1*sin(x) then the measure should be higher than the> above.> Could anyone here help me to find a solution to this?> Anon.If x(t) is the value of time series at time t, you could define thesmoothness of the time series between times t1 and t2 as|x(t2) - x(t1)| / (sum(i=t1+1,t2)(|x(t) - x(t-1)|))In words, this is the ratio of (total displacement from the beginningto the end of the time period) to (sum of the absolute one-periodchanges). The ratio is bounded by the values 0 (jagged) and 1(smooth).I think technical analysts of the financial markets have used thismeasure to determine if a market is trending. === optimization problem, and would like to find some simple intros,tutorials, coding examples, or someplace to ask basic questions. Anewsgroup would be great, but I doubt there is one just on ampl. I'vefound a number of code examples on the web, but they're too complicatedto serve as tutorials.(An example of what I need to start with: I will have a fairly largematrix in the data set -- too large to enter the data explicitly. Butit is sparse, and it would be easy to write something like matlab codeto loop thru the rows/columns and fill in the non-zero data. But allthe ampl code examples I can find just have the data entered manually,not by code.)Ben Crain === Subject: Re: Help on interval arithmetics package for lisp> I'd be grateful for pointers as to where I could find an> interval-arithmetic package for lisp.Nobody has responded to you yet. After all these years, is yourwhen doing a Google search. (Google didn't exist when you posted it!!) === Subject: Re: How to solve recurrence relation?>public static int power2(int base, int exponent) {> if (exponent == 0)> return 1;> if (exponent == 1)> return base;> int part = power2(base, exponent/2);> if (exponent % 2 == 0)> return part * part;> return base * part * part;>} The basic idea of most recursive programming is you handle the simplecases, and you simplify any complex case a tad and feed it back intothe same process. As long as you do a little simplification at eachstage, eventually the process will unwind when it hits one of thetrivial cases.In your example, the simplification is cutting the exponent in halfeach go round. This eventually gets to 1. There are many copies ofpart -- one for each incarnation of power2.Try tracing a simple case, see http://mindprod.com/jgloss/trace.htmlor just doing it on paper.Canadian Mind Products, Roedy Green.Coaching, problem solving, economical contract programming. See http://mindprod.com/jgloss/jgloss.html for The Java Glossary. === Subject: Gambling QuestionHypothetical:You're a gambler, and a good one. You bet professional baseball, butonly on home underdogs, and primarily low odds ones at that (avg.+1.10 or paying 110% of whatever value you risk on any given game).Whereas these particular underdogs have only avgd. winning about 45%of the time over the last several years, your picks win at about a 70%rate. Your system dictates that you place bets on avg. only 2.5 to 3days per week. Mostly you bet only one game on those days where you dobet, but occasionally you will bet 2 or even 3 games on a given day.Even so, you bet a total of slightly less than 4 games per week.Up till now, you have been evaluating the new MLB (Major LeagueBaseball) season without placing any bets, but now you have all theinformation necessary, and are ready to begin betting. You plan toquit betting at the All-Star break in mid-July, giving youapproximately 12 weeks in which to make money. (The second half of theseason doesn't seem to conform to your system.) Your beginning stakeis $3000, and you plan to increase or decrease the size of your betsdepending on the total amt. of money you have won or lost at the timeof any given bet.What is the optimum percentage of your stake (the stake you have atthe time of the bet) which you should bet on those days you have only1 game, days where you have 2 games, and days where you have 3 gamesto bet? (This assumes you will be betting less per game on those dayswhere you have multiple games on the menu.)What would be the anticipated winnings for the season predicated onthe above? === Subject: Re: Gambling QuestionQuestion; You win 70 % of the time, and get 110% of what is bet? Then the30% of the time you loose everything a 0% of what is bet. So the weightedaverage return overall is 77%, or you are loosing 23% of your money no materhow you structure your betting?Sounds reasonable.> Hypothetical:> You're a gambler, and a good one. You bet professional baseball, but> only on home underdogs, and primarily low odds ones at that (avg.> +1.10 or paying 110% of whatever value you risk on any given game).> Whereas these particular underdogs have only avgd. winning about 45%> of the time over the last several years, your picks win at about a 70%> rate. Your system dictates that you place bets on avg. only 2.5 to 3> days per week. Mostly you bet only one game on those days where you do> bet, but occasionally you will bet 2 or even 3 games on a given day.> Even so, you bet a total of slightly less than 4 games per week.> Up till now, you have been evaluating the new MLB (Major League> Baseball) season without placing any bets, but now you have all the> information necessary, and are ready to begin betting. You plan to> quit betting at the All-Star break in mid-July, giving you> approximately 12 weeks in which to make money. (The second half of the> season doesn't seem to conform to your system.) Your beginning stake> is $3000, and you plan to increase or decrease the size of your bets> depending on the total amt. of money you have won or lost at the time> of any given bet.> What is the optimum percentage of your stake (the stake you have at> the time of the bet) which you should bet on those days you have only> 1 game, days where you have 2 games, and days where you have 3 games> to bet? (This assumes you will be betting less per game on those days> where you have multiple games on the menu.)> What would be the anticipated winnings for the season predicated on> the above? === Subject: Re: DEBUGing help needed in FORTRAN coding>Hi there,>I am using a Finite Element Method FORTRAN package in my computation.>With certain geometry (square with small hole), when I use a mesh with>node number less than 400, it works very well, but when I increase the>node number to 1000, it blows up. >Complier said stack over flow. Are you sure that the compiler, rather than the OS, is the one whosaid that ?Did you try limit stack unlimited ?>My question is: how can I fix this kind of bug in FORTRAN? Like how>can I find which stack is full and how to make it larger?>Toby === Subject: Matlab and waveletsI'm trying to determine if it is possible to employ quadratic spline wavelets in Matlab. Can anyone help?Gib === Subject: Re: c code for low pass FIR filter by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i3JAXhT17257;>i need source code for a digital low pass filter in C for>implementation on SHARC EZ-KIT.i need soft copy of filter code(c lang).please send me.thank u, === Subject: Re: Question: to solve nonlinear ODE by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i3JAXhf17231;Please help me solve:L''+aL'-b/L=0This is a nonlinear initial value problem.Initial conditions:L'(0)=0; L(0)=L0L is a function of time.I need a code of Matlab. === Subject: Re: Question: to solve nonlinear ODE >Please help me solve: >L''+aL'-b/L=0 >This is a nonlinear initial value problem. >Initial conditions: >L'(0)=0; L(0)=L0 >L is a function of time. >I need a code of Matlab.translate to a first oder system using y(1)=L and y(2)=L' then use ode45 for example. === Subject: Re: Question: to solve nonlinear ODE> Please help me solve:> L''+aL'-b/L=0> This is a nonlinear initial value problem.> Initial conditions:> L'(0)=0; L(0)=L0> L is a function of time.> I need a code of Matlab.Hi !Maple V R4 solves this problem, but without ICs. To complete the 2 solutionsreturned one has to be able to solve the following two equations respectivelyfor _C1 (a constant of integration) :-1/2*(-2*b*ln(u)+2*_C1)^(1/2)*u*Pi^(1/2)*2^(1/2)*a*(erf(1/2*2 ^(1/2)/b^(1/2)*(-2*b*ln(u)+2*_C1)^(1/2))-erf(1/2*2^(1/2)/b^(1/ 2)*(-2*b*ln(L0)+2*_C1)^(1/2)))/b^(1/2)*exp(-(b*ln(u)-_C1)/b)-( -2*b*ln(u)+2*_C1)^(1/2)=L0*a;-1/2*(-2*b*ln(u)+2*_C1)^(1/2)*u* Pi^(1/2)*2^(1/2)*a*(erf(1/2*2^(1/2)/b^(1/2)*(-2*b*ln(u)+2*_C1) ^(1/2))-erf(1/2*2^(1/2)/b^(1/2)*(-2*b*ln(L0)+2*_C1)^(1/2)))/b^ (1/2)*exp(-(b*ln(u)-_C1)/b)+(-2*b*ln(u)+2*_C1)^(1/2)=L0*a; / 1/2 1/2 1/2 |- 1/2 (-2 b ln(u) + 2 _C1) u Pi 2 a | | 1/2 1/2 2 (-2 b ln(u) + 2 _C1) erf(1/2 ----------------------------) 1/2 b 1/2 1/2 2 (-2 b ln(L0) + 2 _C1) | b ln(u) - _C1 - erf(1/2 -----------------------------)| exp(- -------------) 1/2 | b b / / 1/2 1/2 / b - (-2 b ln(u) + 2 _C1) = L0 a / / 1/2 1/2 1/2 |- 1/2 (-2 b ln(u) + 2 _C1) u Pi 2 a | | 1/2 1/2 2 (-2 b ln(u) + 2 _C1) erf(1/2 ----------------------------) 1/2 b 1/2 1/2 2 (-2 b ln(L0) + 2 _C1) | b ln(u) - _C1 - erf(1/2 -----------------------------)| exp(- -------------) 1/2 | b b / / 1/2 1/2 / b + (-2 b ln(u) + 2 _C1) = L0 a /Chris === Subject: function generationDo anyone know of a program that can generate functionsbased on specified input?Ex:Given f(a1)=b1, f(a2)=b2, f(a3)=b3Generate the general function(s) f().Geir Atle, http://heim.ifi.uio.no/gahegsvo/ If life had an answer, I wouldn't even know the question. === Subject: Re: function generation> Do anyone know of a program that can generate functions> based on specified input?> Ex:> Given f(a1)=b1, f(a2)=b2, f(a3)=b3> Generate the general function(s) f().> --> Geir Atle, http://heim.ifi.uio.no/gahegsvo/> If life had an answer, I wouldn't even know the question.Please see www.estlab.commail Arto.Huttunen@estlabdotcom === Subject: Re: function generation> Do anyone know of a program that can generate functions> based on specified input?> Ex:> Given f(a1)=b1, f(a2)=b2, f(a3)=b3> Generate the general function(s) f().> --> Geir Atle, http://heim.ifi.uio.no/gahegsvo/> If life had an answer, I wouldn't even know the question.Not sure I am understanding your question. What you probably want is an interpolating polynomial/function. given a vector of points and function values at those points, how to get the function value at another point. Is that correct?For the problem you specified, I could think of a polynomial as followsf(x) = b1 * (x-a2) * (x-a3) / (a1 - a2) * (a1 - a3) + b2 * (x-a1) * (x-a3) / (a2 - a1) * (a2 - a3) + b3 * (x-a1) * (x-a2) / (a3 - a1) * (a3 - a2)Standard methods exist, if you want the polynomial in say sin and cosines etc...Home page: http://www.people.cornell.edu/pages/kk288/index.htmlFluid's group: http://groups.yahoo.com/group/flumech/ === Subject: Re: function generation Do anyone know of a program that can generate functions> based on specified input?> Ex:> Given f(a1)=b1, f(a2)=b2, f(a3)=b3> Generate the general function(s) f().> --> Geir Atle, http://heim.ifi.uio.no/gahegsvo/> If life had an answer, I wouldn't even know the question.> Not sure I am understanding your question. What you probably want is an> interpolating polynomial/function. given a vector of points and function> values at those points, how to get the function value at another point.> Is that correct?> For the problem you specified, I could think of a polynomial as follows> f(x) = b1 * (x-a2) * (x-a3) / (a1 - a2) * (a1 - a3)> + b2 * (x-a1) * (x-a3) / (a2 - a1) * (a2 - a3)> + b3 * (x-a1) * (x-a2) / (a3 - a1) * (a3 - a2)> Standard methods exist, if you want the polynomial in say sin and> cosines etc...> Home page: http://www.people.cornell.edu/pages/kk288/index.html> Fluid's group: http://groups.yahoo.com/group/flumech/What I want is to find a function based on known input to and outputfrom the function.Say we have f(1)=1, f(2)=3, f(3)=6, f(4)=10, ...One fuction to match this input and output would be:f(x) = x * (x+1) / 2Now, say we didn't know f(x) was the correct function.Is there a program we could feed the input and output to sothat it would generate f(x) for us? === Subject: Re: function generation> What I want is to find a function based on known input to and output> from the function.> Say we have f(1)=1, f(2)=3, f(3)=6, f(4)=10, ...> One fuction to match this input and output would be:> f(x) = x * (x+1) / 2> Now, say we didn't know f(x) was the correct function.The problem is that, based on that data, you *still* don't know that f(x)is the correct function. The correct function could also bef(x) = x*(x+1)/2 + 1000000*(x-1)*(x-2)*(x-3)*(x-4)> Is there a program we could feed the input and output to so> that it would generate f(x) for us?Yes, assuming that both of the above choices of f are equally good inyour eyes. If they aren't equally good, then you'll need to add somethingto your problem statement to discriminate between good and bad choices off.--- === Subject: Re: function generation What I want is to find a function based on known input to and output> from the function.> Say we have f(1)=1, f(2)=3, f(3)=6, f(4)=10, ...> One fuction to match this input and output would be:> f(x) = x * (x+1) / 2> Now, say we didn't know f(x) was the correct function.> The problem is that, based on that data, you *still* don't know that f(x)> is the correct function. The correct function could also be> f(x) = x*(x+1)/2 + 1000000*(x-1)*(x-2)*(x-3)*(x-4)> Is there a program we could feed the input and output to so> that it would generate f(x) for us?> Yes, assuming that both of the above choices of f are equally good in> your eyes. If they aren't equally good, then you'll need to add something> to your problem statement to discriminate between good and bad choices of> f.> ---> Roy StognerYes, the above choices are equally good. What program would do the trick?Geir Atle, http://heim.ifi.uio.no/gahegsvo/ If life had an answer, I wouldn't even know the question. === Subject: Re: function generation>What I want is to find a function based on known input to and output>from the function.>Say we have f(1)=1, f(2)=3, f(3)=6, f(4)=10, ...>One fuction to match this input and output would be:>f(x) = x * (x+1) / 2>Now, say we didn't know f(x) was the correct function.> The problem is that, based on that data, you *still* don't know that f(x)> is the correct function. The correct function could also be> f(x) = x*(x+1)/2 + 1000000*(x-1)*(x-2)*(x-3)*(x-4)>Is there a program we could feed the input and output to so>that it would generate f(x) for us?> Yes, assuming that both of the above choices of f are equally good in> your eyes. If they aren't equally good, then you'll need to add something> to your problem statement to discriminate between good and bad choices of> f.> ---> Roy Stogner> Yes, the above choices are equally good. What program would do the trick?If canned software is OK, try InterpolatingFunction[domain,table] in Mathematica. Are you by chance related to Kjell Magne Mathisen at Trondheim? === Subject: Re: function generation> What I want is to find a function based on known input to and output> from the function.> Say we have f(1)=1, f(2)=3, f(3)=6, f(4)=10, ...> One fuction to match this input and output would be:> f(x) = x * (x+1) / 2> Now, say we didn't know f(x) was the correct function.> The problem is that, based on that data, you *still* don't know that f(x)> is the correct function. The correct function could also be> f(x) = x*(x+1)/2 + 1000000*(x-1)*(x-2)*(x-3)*(x-4)> Is there a program we could feed the input and output to so> that it would generate f(x) for us?> Yes, assuming that both of the above choices of f are equally good in> your eyes. If they aren't equally good, then you'll need to add something> to your problem statement to discriminate between good and bad choices of> f.> Yes, the above choices are equally good.Wow, I admit I wasn't expecting that answer. ;-) Usually people want tofit curves to four points so they can extrapolate a fifth, and becomeconcerned when they discover that the simplest extrapolation technique cango wild.Surely you must have some additional requirements for f; otherwise why notjust let your function f be a piecewise constant around each data point?What application is this for?> What program would do the trick?I'm all in favor of reusing code, but this is simple enough that you'llwant to rewrite it yourself. If you have (or can quickly write) softwarefor multiplying polynomial objects, then one way to do it is to generalizethe formula gave.If you really don't have any stricter requirements on f(), why not justmake your function be piecewise linear? If your data is f(a_i) = b_i,then let f(x) := b_i + (b_(i+1) - b_i) * (a_(i+1) - a_i) for a_i < x < a_(i+1),f(x) := b_1 for x < a_1,f(x) := b_n for x > a_n--- === Subject: Re: function generation>Do anyone know of a program that can generate functions>based on specified input?>Ex:>Given f(a1)=b1, f(a2)=b2, f(a3)=b3>Generate the general function(s) f().>-->Geir Atle, http://heim.ifi.uio.no/gahegsvo/> If life had an answer, I wouldn't even know the question.>Not sure I am understanding your question. What you probably want is an>interpolating polynomial/function. given a vector of points and function>values at those points, how to get the function value at another point.>Is that correct?>For the problem you specified, I could think of a polynomial as follows>f(x) = b1 * (x-a2) * (x-a3) / (a1 - a2) * (a1 - a3)> + b2 * (x-a1) * (x-a3) / (a2 - a1) * (a2 - a3)> + b3 * (x-a1) * (x-a2) / (a3 - a1) * (a3 - a2)>Standard methods exist, if you want the polynomial in say sin and>cosines etc...>Home page: http://www.people.cornell.edu/pages/kk288/index.html>Fluid's group: http://groups.yahoo.com/group/flumech/> What I want is to find a function based on known input to and output> from the function.> Say we have f(1)=1, f(2)=3, f(3)=6, f(4)=10, ...> One fuction to match this input and output would be:> f(x) = x * (x+1) / 2> Now, say we didn't know f(x) was the correct function.> Is there a program we could feed the input and output to so> that it would generate f(x) for us?> --> Geir Atle, http://heim.ifi.uio.no/gahegsvo/> If life had an answer, I wouldn't even know the question.(a) for any finite number of (x, y) pairs that you have to build a function out of, there's an infinite number of functions that you can build that will map all of the x's to the correct y.(b) for any function you choose to to map x onto y that works for the finite set, I can find an x, y pair that doesn't fit.There's a plethora of ways to make approximations to the correct function for a given set of (x, y) pairs, ranging from least-squares fits to some assumed function to neural nets, but all of them require some starting points supplied by a knowledgeable human being.http://www.wescottdesign.com === Subject: Curve Approximation by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i3JE0rp15227;My problem is to write a program in C language that approximates apolynomial(up to degree 7) curve from a thermocoupleboard(x:temperature; y:voltage).I'm searching for a detailed algorythm or a C source that does thejob.Guilaume Le GoffStudent in electronics and electrotechnicsat ESIEE(french engineers school) === Subject: Re: Curve ApproximationThe least squares fit for a 7th degree polynomial can be extremely ill-conditioned. You should use something like the algorithms at http://lib.stat.cmu.edu/apstat/274 http://lib.stat.cmu.edu/apstat/75to avoid computing X'XJerry> My problem is to write a program in C language that approximates a> polynomial(up to degree 7) curve from a thermocouple> board(x:temperature; y:voltage).> I'm searching for a detailed algorythm or a C source that does the> job. === Subject: Re: Curve Approximation> My problem is to write a program in C language that approximates a> polynomial(up to degree 7) curve from a thermocouple> board(x:temperature; y:voltage).> I'm searching for a detailed algorythm or a C source that does the> job.> Guilaume Le Goff> Student in electronics and electrotechnics> at ESIEE(french engineers school)Here are some well meaning linesIn some cases 7th degree might be reasonable.In regression analysis x,x^2,x^3... are strongly correlated and causenumerical problems in calculation.Take care of some precautions.Use normalized variables for all x,x^2,x^3...Use double precision in calculationsUse conjugate gradient algorithm to solve the coefficientsMoreUse large set of x,y pairs much more than 7Make experiments with 3, 5, 7 degreesand see if you can have a better approximation.Take some of the x,y pairs away and see if thecoefficients remain essentially the same.If you have coefficient values like 7.5E16 for x^4 and -1.234E17 for x^5then you are apparently modeling the signal noise.Arto.Huttunen@estlabdotcom === Subject: Re: Curve Approximation> My problem is to write a program in C language > that approximates a polynomial(up to degree 7) curve> from a thermocouple board(x:temperature; y:voltage).> I'm searching for a detailed algorithm or a C source > that does the job.Take a look atThe C++ Scalar, Vector, Matrix and Tensor [class] Library: http://www.netwood.net/~edwin/svmtl/ > cat svmt/examples/polynomial.cc/*The C++ Scalar, Vector, Matrix and Tensor classes.Copyright (C) 1998 E. Robert TisdaleThis file is part of The C++ Scalar, Vector, Matrix and Tensor classes.This library is free software which you can redistribute and/or modifyunder the terms of the GNU Library General Public Licenseas published by the Free Software Foundation;either version 2, or (at your option) any later version.This library is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY; without even the implied warrantyof MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.See the GNU General Public License for more details.You should have received a copy of the GNU Library General Public Licensealong with this library. If not, write to the Free Software Foundation,Inc., 675 Mass Ave, Cambridge, MA 02139, USA.Written by E. Robert Tisdale*//*This little program fits a polynomial y(x) = w_0 + w_1*x + w_2*x^2 + ... + w_n*x^nof order n to a set of m > n data pairs {x_i, y_i}from standard input and prints the coefficients to standard output. $ make polynomial $ polynomial 2 < data w = -3.44668 1.01307 0.499563 $ gnuplot gnuplot> y(x) = -3.44668 + 1.01307*x + 0.499563*x*x gnuplot> plot data, y(x)*/#include x && cin > y) { // {x_i, y_i} for (Offset j = 0; j < n; ++j) u[j+1] = x*u[j]; // u_ij = (x_i)^j S += u.submatrix().t().dot(); // S += ((u_i)^T)u_i v += y*u; // v += y_i*u_i } doubleVector w = solve(v, S); // solve v = wS for w cout << w =n << w; return 0; } === Subject: Re: Curve Approximation >My problem is to write a program in C language that approximates a >polynomial(up to degree 7) curve from a thermocouple >board(x:temperature; y:voltage). >I'm searching for a detailed algorythm or a C source that does the >job. >Guilaume Le Goff >Student in electronics and electrotechnics >at ESIEE(french engineers school) look here:http://plato.la.asu.edu/topics/problems/nlolsq.htmlthere is a list of codes for doing least squares, also in c.but as mentioned by someone else already, using a polynomialof order 7 in this application is more than questionable. first look at the shape of your data. maybe some exponential fit is much moreadequate here. === Subject: Re: Curve Approximation> My problem is to write a program in C language that approximates a> polynomial(up to degree 7) curve from a thermocouple> board(x:temperature; y:voltage).> I'm searching for a detailed algorythm or a C source that does the> job. I would be deeply suspicious of _anything_ that uses a 7th-order polynomial to fit to established data. Polynomial curve fits tend to be very badly behaved when you run off the end of the calibrated portion. It's much better to start with a model that's close to what's going on physically, fit to that, then use it. When you run off the end of your calibrated data you'll be doing much more intelligent extrapolation.http://www.wescottdesign.com === Subject: System of ODEsMy aim is the solution of a huge, sparse system of ordinary firstorder differential equations of the general formdy/dt = A y = (B + c D) yHere, A and B are Hermitian matrices and D is a diagonal matrix. Thesolution is to be obtained for a series of equally spaced values ofthe real parameter c. How can I make use of the special structure ofthe problem? Assuming a solution via the matrix exponent is attempted,is it possible to use the eigenvectors obtained for one value of c toaid the calculation of y for subsequent values of c? By the way, cspans several orders of magnitude and, hence, c D may not always beregarded a small perturbation. However, adjacent values of c differonly slightly.Any ideas?Johan === Subject: Re: System of ODEs> My aim is the solution of a huge, sparse system of ordinary first> order differential equations of the general form> dy/dt = A y = (B + c D) y> Here, A and B are Hermitian matrices and D is a diagonal matrix. The> solution is to be obtained for a series of equally spaced values of> the real parameter c. How can I make use of the special structure of> the problem? Assuming a solution via the matrix exponent is attempted,> is it possible to use the eigenvectors obtained for one value of c to> aid the calculation of y for subsequent values of c? By the way, c> spans several orders of magnitude and, hence, c D may not always be> regarded a small perturbation. However, adjacent values of c differ> only slightly.> Any ideas?I assume D is a positive matrix (not really necessary, but helpful), andthat B and D are time-independent. Then Step 1: transform the dependent variable, x = sqrt(D) y; Step 2: now the equation is dx/dt = Cx + cIx where I = unit matrix and C = [sqrt(D)]^{-1} B [sqrt(D)]^{-1}; Step 3: define z = x exp(-ct) the equation is now dz/dt = Cz Step 4: solve by inspection in terms of exponentiated matrix C (C is still Hermitian).Julian V. NobleProfessor Emeritus of Physicsjvn@lessspamformother.virginia.edu ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ For there was never yet philosopher that could endure the toothache patiently. -- Wm. Shakespeare, Much Ado about Nothing. Act v. Sc. 1. === Subject: Re: System of ODEs> My aim is the solution of a huge, sparse system of ordinary first> order differential equations of the general form> dy/dt = A y = (B + c D) y> Here, A and B are Hermitian matrices and D is a diagonal matrix. The> solution is to be obtained for a series of equally spaced values of> the real parameter c. How can I make use of the special structure of> the problem? Assuming a solution via the matrix exponent is attempted,> is it possible to use the eigenvectors obtained for one value of c to> aid the calculation of y for subsequent values of c? By the way, c> spans several orders of magnitude and, hence, c D may not always be> regarded a small perturbation. However, adjacent values of c differ> only slightly.> Any ideas?> I assume D is a positive matrix (not really necessary, but helpful), and> that B and D are time-independent. Then> Step 1: transform the dependent variable, x = sqrt(D) y;Unfortunately D is not in general positive, though, B and D are indeedtime-independent as you assumed. Even worse, D is not in generalinvertible. > Step 2: now the equation is dx/dt = Cx + cIx where I = unit matrix> and C = [sqrt(D)]^{-1} B [sqrt(D)]^{-1};Sorry for my ignorance, but I simply can not comprehend your result.I got dx/dt = (D^{1/2} B D^{-1/2} + c D) x, but maybe I need to getsome more coffee first. Please clarify.> Step 3: define z = x exp(-ct) the equation is now dz/dt = Cz> Step 4: solve by inspection in terms of exponentiated matrix C> (C is still Hermitian).Johan === Subject: Re: System of ODEs>My aim is the solution of a huge, sparse system of ordinary first>order differential equations of the general form>dy/dt = A y = (B + c D) y>Here, A and B are Hermitian matrices and D is a diagonal matrix. The>solution is to be obtained for a series of equally spaced values of>the real parameter c. How can I make use of the special structure of>the problem? Assuming a solution via the matrix exponent is attempted,>is it possible to use the eigenvectors obtained for one value of c to>aid the calculation of y for subsequent values of c? By the way, c>spans several orders of magnitude and, hence, c D may not always be>regarded a small perturbation. However, adjacent values of c differ>only slightly.>Any ideas?> I assume D is a positive matrix (not really necessary, but helpful), and> that B and D are time-independent. Then> Step 1: transform the dependent variable, x = sqrt(D) y;> Unfortunately D is not in general positive, though, B and D are indeed> time-independent as you assumed. Even worse, D is not in general> invertible.> Step 2: now the equation is dx/dt = Cx + cIx where I = unit matrix> and C = [sqrt(D)]^{-1} B [sqrt(D)]^{-1};> Sorry for my ignorance, but I simply can not comprehend your result.> I got dx/dt = (D^{1/2} B D^{-1/2} + c D) x, but maybe I need to get> some more coffee first. Please clarify.> Step 3: define z = x exp(-ct) the equation is now dz/dt = Cz> Step 4: solve by inspection in terms of exponentiated matrix C> (C is still Hermitian).Not your ignorance, my bad algebra. Don't know what I could have beenthinking of (probably a bad case of mental languish), since this is justtime-dependent perturbation theory in quantum mechanics.What I meant was y = exp[cDt] x dx/dt = exp[-cDt] B exp[cDt] x = C(t) x .This may or may not be helpful, since the solution is the time-orderedproduct x(t) = T exp{ int_0^t {ds,C(s)} } x(0) ,which is generally a bear to compute.If the commutator of B and D is not a constant (times the unit matrix),you also can't use the Baker-Campbell-Hausdorff theorem, exp(A+B) = exp(A) exp([B,A]/2) exp(B) ,where [B,A] = BA - AB . If any higher commutator vanishesthere is a generalization of the BCH theorem that might help.The generalization was discussed in The Journal of MathematicalPhysics some years ago. Personally I have only found the simplestversion of BCH helpful.Sorry, that's the best I can do. Julian V. NobleProfessor Emeritus of Physicsjvn@lessspamformother.virginia.edu ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ For there was never yet philosopher that could endure the toothache patiently. -- Wm. Shakespeare, Much Ado about Nothing. Act v. Sc. 1. === Subject: Re: System of ODEs> My aim is the solution of a huge, sparse system of ordinary first> order differential equations of the general form> dy/dt = A y = (B + c D) y> Here, A and B are Hermitian matrices and D is a diagonal matrix. The> solution is to be obtained for a series of equally spaced values of> the real parameter c. How can I make use of the special structure of> the problem? Assuming a solution via the matrix exponent is attempted,> is it possible to use the eigenvectors obtained for one value of c to> aid the calculation of y for subsequent values of c? By the way, c> spans several orders of magnitude and, hence, c D may not always be> regarded a small perturbation. However, adjacent values of c differ> only slightly.> Any ideas?> I assume D is a positive matrix (not really necessary, but helpful), and> that B and D are time-independent. Then> Step 1: transform the dependent variable, x = sqrt(D) y;> Step 2: now the equation is dx/dt = Cx + cIx where I = unit matrix> and C = [sqrt(D)]^{-1} B [sqrt(D)]^{-1};and since [C,cI]=0, x(t)=exp(ct)exp(Ct)x(0). If [B,D]=0 theny(t)=exp(cDt)exp(Bt). So c is not the problem but rather having tocalculate exp(Ct) or exp(Bt) when A is huge and sparse (selfadjointness isof no advantage except that the eigenvalues of A are real) is a problem andis usually done via Kyrov subspace techniques.E&OECiao,Gerry T. === Subject: Re: System of ODEs>My aim is the solution of a huge, sparse system of ordinary first>order differential equations of the general form>dy/dt = A y = (B + c D) y>Here, A and B are Hermitian matrices and D is a diagonal matrix. The>solution is to be obtained for a series of equally spaced values of>the real parameter c. How can I make use of the special structure of>the problem? Assuming a solution via the matrix exponent is attempted,>is it possible to use the eigenvectors obtained for one value of c to>aid the calculation of y for subsequent values of c? By the way, c>spans several orders of magnitude and, hence, c D may not always be>regarded a small perturbation. However, adjacent values of c differ>only slightly.>Any ideas?> I assume D is a positive matrix (not really necessary, but helpful), and> that B and D are time-independent. Then> Step 1: transform the dependent variable, x = sqrt(D) y;> Step 2: now the equation is dx/dt = Cx + cIx where I = unit matrix> and C = [sqrt(D)]^{-1} B [sqrt(D)]^{-1};> and since [C,cI]=0, x(t)=exp(ct)exp(Ct)x(0). If [B,D]=0 then> y(t)=exp(cDt)exp(Bt). So c is not the problem but rather having to> calculate exp(Ct) or exp(Bt) when A is huge and sparse (selfadjointness is> of no advantage except that the eigenvalues of A are real) is a problem and> is usually done via Kyrov subspace techniques.Indeed, B and D to not commute. I forgot to mention that before.consider using expokit in order to calculate exp(A t)*y_0 before.Still I hoped I could somehow make use of previous results whenstepping c.Johan === Subject: Re: System of ODEs>My aim is the solution of a huge, sparse system of ordinary first>order differential equations of the general form>dy/dt = A y = (B + c D) y>Here, A and B are Hermitian matrices and D is a diagonal matrix.The>solution is to be obtained for a series of equally spaced values of>the real parameter c. How can I make use of the special structureof>the problem? Assuming a solution via the matrix exponent isattempted,>is it possible to use the eigenvectors obtained for one value of cto>aid the calculation of y for subsequent values of c? By the way, c>spans several orders of magnitude and, hence, c D may not always be>regarded a small perturbation. However, adjacent values of c differ>only slightly.>Any ideas?>I assume D is a positive matrix (not really necessary, but helpful),and>that B and D are time-independent. Then> Step 1: transform the dependent variable, x = sqrt(D) y;> Step 2: now the equation is dx/dt = Cx + cIx where I = unit matrix> and C = [sqrt(D)]^{-1} B [sqrt(D)]^{-1};> and since [C,cI]=0, x(t)=exp(ct)exp(Ct)x(0). If [B,D]=0 then> y(t)=exp(cDt)exp(Bt). So c is not the problem but rather having to> calculate exp(Ct) or exp(Bt) when A is huge and sparse (selfadjointnessis> of no advantage except that the eigenvalues of A are real) is a problemand> is usually done via Kyrov subspace techniques.> Indeed, B and D to not commute. I forgot to mention that before.> consider using expokit in order to calculate exp(A t)*y_0 before.> Still I hoped I could somehow make use of previous results when> stepping c.the vector exp(Ct)D^1/2y(0) (computationally expensive) and haveD^-1/2exp(ct) (computationally inexpensive) operate on it over the range ofc. Doesn't come much better than that, :-).E&OECiao,Gerry T. === Subject: Re: System of ODEs> My aim is the solution of a huge, sparse system of ordinary first> order differential equations of the general form> dy/dt = A y = (B + c D) y> Here, A and B are Hermitian matrices and D is a diagonal matrix. The> solution is to be obtained for a series of equally spaced values of> the real parameter c. How can I make use of the special structure of> the problem? Assuming a solution via the matrix exponent is attempted,> is it possible to use the eigenvectors obtained for one value of c to> aid the calculation of y for subsequent values of c? By the way, c> spans several orders of magnitude and, hence, c D may not always be> regarded a small perturbation. However, adjacent values of c differ> only slightly.> Any ideas?> I assume D is a positive matrix (not really necessary, but helpful), and> that B and D are time-independent. Then> Step 1: transform the dependent variable, x = sqrt(D) y;> Step 2: now the equation is dx/dt = Cx + cIx where I = unit matrix> and C = [sqrt(D)]^{-1} B [sqrt(D)]^{-1}; ^ sorry, forgot here> Step 3: define z = x exp(-ct) the equation is now dz/dt = Cz> Step 4: solve by inspection in terms of exponentiated matrix Ct ^ i.e. you want exp[Ct]> (C is still Hermitian). Sorry for the typos.Julian V. NobleProfessor Emeritus of Physicsjvn@lessspamformother.virginia.edu ^^^^^^^^^^^^^^^^^^http://galileo.phys.virginia.edu/~jvn/ For there was never yet philosopher that could endure the toothache patiently. -- Wm. Shakespeare, Much Ado about Nothing. Act v. Sc. 1. === Subject: scaat kalman filter with simultaneous observationsI am trying to implement a scaat kalman filter,and will get several measurements from a cameraall at the same frame time. I can see how toimplement this with a different time for eachmeasurement (even when they come out of sequence),but the Jacobian seems impossible to calculatewhen the denominator is zero. The delta ofthe state matrix elements is zero.Is is correct to calculate the Jacobian numerically?How can I do this when the measurement is atthe same time?See http://www.cs.unc.edu/~welch/media/pdf/scaat_ dissertation.pdfchapter 4 equation 4.5 and 4.12Jack === Subject: Re: sparse matrix eigenvalue problem > (see below) > aaah, you want a general sparse QR iteration? I think this does not exist because > the extreme fill in will destroy a general sparsitiy pattern. these exist > codes for qr-iteration for banded matrices and you might try to first transform > _once_ your matrix into a banded one with a reasonable wide band. >This works only in the symmetric case, and then one can directly reduce >to tridiagonal form. no, unfortunately. if you want to achieve directly' the tridiagonal band form, then you get an complete fill in intermediately, absolutely disastrous. this is the reason why h.r.schwarz invented his quite tricky chasing for the bandwidth-reduction by unitary similarity transforms. also this is costly peter >In general, QR only preserves the Hessenberg form, and doubles >the 'other' bandwith in every step, unless this is prevented by symmetry. >Thus doing the band reduction is not useful in the nonsymmetric case. >But one can reduce and then use Laguerre iteration or inverse iteration >with Rayleigh quotient shift ... >Arnold Neumaier === Subject: Re: Best fit of rectangles in a bigger rectangle already having prefilled rectangles >i need a alogrithm for computing the the locations of the rectangles >that can be best fitted in another rectangle. The other rectangle is >already been filled with rectangles with permanent non-changeable >locations and size. >The Algorithm should be able to fit other rectangles given their size >and to fit as close as to the other rectangles if possible without >colliding with other each other. >any programs or logic would be appreciates. please let me know if u >need more details. >rajaonce again the 2D cutting stock problem. yes, there exist algorithms forapplication (designing a webpage), this is a big overshoot. look on the screen and do it yourself === Subject: Re: SVD performance, Matlab vs. Lapack > I'm making a lot (many thousands) of singular value decompositions >of > medium-sized matrices ( e.g. 100 x 50 ) >What do you need these for? >I want to experiment with regularizing these matrices before >inversion. The physical problem relates to acoustics. Anyway, the >subject here is performance, how to efficiently decompose complex >matrices. since MATLAB uses LAPACK, if you see _your_ lapck slower than matlab menas that you did no use optimization plus an optimized BLAS(since LAPACK rleis on an efficient BLAS). clearly MATLAB is shipped with anoptimized BLAS. If you miss one on your system, download ATLAS from netliband install it yourself. It does a very good job.and of course, use full optimization in your compiler call === Subject: Re: SVD performance, Matlab vs. Lapack> I guess I'll have to dive into the vast world of fortran compilation. > Sadly enough, since I don't expect to improve performance considerably> that way. But who knows? I had hoped for some SVD expert to say:Oh, didn't you know that you can make an adequate performance-cheap> SVD in this-and-this way?. I tried to use less precision (cgesvd> instead of zgesvd) in order to minimize iterative steps, but> performance was exactly the same.Don't be so glum :). If you use Matlab a lot, I think Fortran 90/95should be a tool in your arsenal. The two languages are not thatdifferent.There is no reason a good Fortran compiler should compute SVD's 4times slower than Matlab, unless you are using very inefficientsoftware (which Lapack with the proper BLAS is not) or Matlab isexploiting some special structure of the matrix that the Fortran codeis not. If Fortran did such linear algebra calculations 4 timesslower, it would be dead. === Subject: Elementary MathematicsWelcome to the All Elementary Mathematics - for people who need help, want to help, or want to find the truth in our onlinediscussions. We can help with all areas of Math. Post your questionshere and someone will respond. www.bymath.com forum www.bymath.com/cgi-bin/ultimatebb.cgi === Subject: choice of program/languageI am looking for the best program/language to write my own developedstatistical procedures in.Up til now everything is programmed in SAS IML.We are looking for another program/language to improve the speed of ourprogram.Other options under consideration are:1. C++2. Matlab3. Mathematica4. ...I have no real experience with any of the above programs/language.Things that are often used: likelihood calculations, matrix calculations(sum, product, kronecker product, inverse,..), optimatization routines (maxlikelihood; possibly constrained), ...Please advice. === Subject: Re: choice of program/language> I am looking for the best program/language to write my own developed> statistical procedures in.> Up til now everything is programmed in SAS IML.> We are looking for another program/language to improve the speed of our> program.> Other options under consideration are:> 1. C++> 2. Matlab> 3. Mathematica> 4. ...> I have no real experience with any of the above programs/language.> Things that are often used: likelihood calculations, matrix calculations> (sum, product, kronecker product, inverse,..), optimatization routines(max> likelihood; possibly constrained), ...> Please advice.Matlab is very good, but very expensive. It has a statistics toolbox (formore money). If money is not a constraint, Matlab is a great choice. Speedmay or may not be an issue depending on your specific needs.IMHO, I would avoid Object Oriented languages (e.g., C++) unless you canidentify clear benefits from using them.There is a huge code base with Fortran (F77, F90 and F95) for yourapplications. Check out netlib.org. === Subject: Re: choice of program/language> I am looking for the best program/language to write my own developed> statistical procedures in.> Up til now everything is programmed in SAS IML.> We are looking for another program/language to improve the speed of our> program.> Other options under consideration are:> 1. C++> 2. Matlab> 3. Mathematica> 4. ...> I have no real experience with any of the above programs/language.> Things that are often used: likelihood calculations, matrix calculations> (sum, product, kronecker product, inverse,..), optimatization routines> (max> likelihood; possibly constrained), ...> Please advice.> Matlab is very good, but very expensive. It has a statistics toolbox (for> more money). If money is not a constraint, Matlab is a great choice. Speed> may or may not be an issue depending on your specific needs.> IMHO, I would avoid Object Oriented languages (e.g., C++) unless you can> identify clear benefits from using them.> There is a huge code base with Fortran (F77, F90 and F95) for your> applications. Check out netlib.org.> JohnI agree with the advice above. For statistics and other topics, aninvaluable source of Fortran 90/95 code is Alan Miller's sitehttp://users.bigpond.net.au/amiller/ . There are links to statisticscodes in Fortran at http://www.dmoz.org/Computers/Programming/Languages/Fortran/ Source_Code/Statistics_and_Econometrics/.A basic problem with C and C++ for numerical work is thatmultidimensional arrays are not built in. Each numerical library seemsto define its own matrix and vector classes, which inhibits reuse. === Subject: Re: choice of program/language> Matlab is very good, but very expensive. It has a statistics toolbox (for> more money). If money is not a constraint, Matlab is a great choice. Speed> may or may not be an issue depending on your specific needs.> IMHO, I would avoid Object Oriented languages (e.g., C++) unless you can> identify clear benefits from using them.> There is a huge code base with Fortran (F77, F90 and F95) for your> applications. Check out netlib.org.> JohnAnother problem with F90 etc., is that there is no free compiler available. Ifc is free, but is not that free as you might want it to be. For example there is a non-commercial unsupported version for linux, What to do if you want to compile the programs on windows/mac environment?I would go for c or c++ where compilers like gcc and g++ are available which are free and portable. The support for this is very good and user community is also large.Aside from that C or Fortran90 is good enough for numerical computations. Stay away from fortran 77, It is one horrible language in the modern world (with all its goto's etc)Home page: http://www.people.cornell.edu/pages/kk288/index.htmlFluid's group: http://groups.yahoo.com/group/flumech/ === Subject: Re: choice of program/language> Another problem with F90 etc., is that there is no free compiler> available...Say, you have a free access to a library, but it only has one book. Downthe road is another library that charges a bit but it has millions ofbooks. Now, that's how you should view Fortran compilers - a bargain atany price!Furthermore, do not get dissuaded by the absence of free F90 compilerssince nearly all F77 compilers already offered (as an extension to thestandard) many advanced features that eventually made it into subsequentstandards. Watcom compiler is but one, and most likely, the finest suchchoice. http://openwatcom.org === Subject: Re: choice of program/language> Aside from that C or Fortran90 is good enough for numerical> computations. Stay away from fortran 77, It is one horrible language in> the modern world (with all its goto's etc)I guess I would have to quibble a bit with this statement. If you want afree compiler, you can get a free F77 compiler.Code in Fortran 77 does not have to use GOTO's. Plenty of F77 code existsthat uses good structured programming practices without a single GOTO.Lapack (which the Matlab computational engine is based on) is an excellentexample of good F77 code. Lapack is probably the best linear algebrasoftware in existence, it's free, and it's written in F77 (although you canuse the Fortran 95 interface if you want).> Home page: http://www.people.cornell.edu/pages/kk288/index.html> Fluid's group: http://groups.yahoo.com/group/flumech/ === Subject: Re: choice of program/languageKris,I personally prefer C++, but you could also choose ADA, becausethey support Object Orientation. If you are prepared to learnhow to design and write OO programs, then it is really worththe effort. On the other hand if you are not interested in OO,and want to write quick programs, that you then throw away,it is not worth the effort to learn and then design in OO.The up side of OO languages is that the programs once writtencan be split up and the parts reused. Also it is easier touse OO libraries from the web eg Blitz++.I am a professional programmer, and do not have extensiveexperience of statistical programs.Hope that helps, Jack> I am looking for the best program/language to write my own developed> statistical procedures in. === Subject: Re: choice of program/language> Kris,> I personally prefer C++, but you could also choose ADA, because> they support Object Orientation. If you are prepared to learn> how to design and write OO programs, then it is really worth> the effort. On the other hand if you are not interested in OO,> and want to write quick programs, that you then throw away,> it is not worth the effort to learn and then design in OO.> The up side of OO languages is that the programs once written> can be split up and the parts reused. Also it is easier to> use OO libraries from the web eg Blitz++. I am looking for the best program/language to write my own developed> statistical procedures in.> Up til now everything is programmed in SAS IML.> We are looking for another program/language to improve the speed of our> program.> Other options under consideration are:> 1. C++> 2. Matlab> 3. Mathematica> 4. ...> I have no real experience with any of the above programs/language.> Things that are often used: likelihood calculations, matrix calculations> (sum, product, kronecker product, inverse,..), optimatization routines (max> likelihood; possibly constrained), ... Don't forget Fortran (you get a personal yecch from me, but it's an option).There are math libraries to do all the matrix stuff you want in C++ or (especially) Fortran, and your file I/O should be faster and less format constrained than with Matlab (and I assume even more so than Mathmatica).In theory using C++ or Fortran should be as fast or faster than Matlab, but they have very carefully optimized routines for the matrix functions; you may not be able to beat them. Anything that you do custom should go at least as fast or faster than what you'd see with Matlab scripting language.Unless you spend a _lot_ of time waiting around for your computer to finish I'd choose the language that you're most comfortable in, or the one you most want to learn.http://www.wescottdesign.com === Subject: Re: choice of program/languageC and/or C++my personal choice of course... === Subject: Re: When should I stop iterative method? by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i3JJbvh01583;>I use iterative methods (Jacobi, Gauss-Seidel, ...) for solvinglinear>equation systems and i want to know when should I stop the iteration?>I know that the condition for stopping has something to do withBanach's>theorem or something like this:>d(x*,x^(k+1))<=q/(1-q)[d(x^(k),x^(k+1)], x* is real solution;>but i dont understand it!>So, what the condition to stop iterative procedure for solving linear>equation systems?In solvingthe system of equations Ax-b = 0,iterarate by whatever method. At each iteration step k evaluate (usingLatex)y_i^{(k)} = mod{ frac{ Sigma a_{ij}x_j^{(k+1)} - b_i}{ Sigmaa_{ij}x_j^{(k)} - b_i}}.If y_i^{(k)} < 1 for each i, as k proceeds, the method is converging.Stop when y_i^{(k)} = 0(or a satisfactory zero tolerance)Note, specifying number of iterations may prematurely terminatethe process; then the wrong conclusion follows - the method does notconverge.O. OtolorinOlabisi Onabanjo UniversityAgo-Iwoye, Nigeria. === Subject: Re: When should I stop iterative method?> At each iteration step k evaluate (using> Latex)> y_i^{(k)} = mod{ frac{ Sigma a_{ij}x_j^{(k+1)} - b_i}{ Sigma> a_{ij}x_j^{(k)} - b_i}}.In other words, componentwise the k+1st residual divided by the kth.> If y_i^{(k)} < 1 for each i, as k proceeds, the method is converging.> Stop when y_i^{(k)} = 0(or a satisfactory zero tolerance)Are you supposing superlinear convergence here? If you have linearconvergence, your y_i^k will converge to a constant between 0 and 1.V.email: lastname at cs utk eduhomepage: cs utk edu tilde lastname === Subject: Logarithmic singularity in integrands by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i3JJc1H01689;I am doing some computations involving computations of integrals withsingular integrands. These integrals appear during application ofboundary integral equations methods to boundary-value problems. Incase of balanced coefficients (Peclet number) well known method basedon splitting integrand function works fine, but in case of singularperturbation (big Peclet number) scheme breaks down, probably, becauseof integrand's behaviour : integrand are delta function - like. Thereare approaches to deal with such problems based on analyticaltransformations rescaling integration region and correctingintegrands in numerical sense. If one has some experience in such problems, please, provide yourfeedback. === Subject: how to select multiple resion of interests in MATLAB roipoly command by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i3K1Fr514132; I have to select multiple polygons in a matlab figure. I think ihave to use roipoly command for that but this command allows me toselect only one polygon at a time and displays it as a new figure. Butin my project i have 50-70 polygons in a single image and i have toselect them all at a time. Any help on this will be greatly appreciated === Subject: Re: Yikes! I've Broken Incommensurability!I dont quite understand some of the words in quotes.Can you make your argument more precise?> 9. This seems, to me, to be quite astounding,> because Incommensurability is nowhere to> be seen.> so, Continuously, without a hint of Incom-> mensurability.> are no longer a Problem - just construce> a Circle with the necessary radious, and, voil?,> there's your Conjugated-pair of Irrationals,> Numbers, but you can address them, Exactly, === Subject: Re: Yikes! I've Broken Incommensurability!> I dont quite understand some of the words in quotes.> Can you make your argument more precise?I will, at a future date.Until then, everyone should remain Skepticalbecause I've not discussed the method Com-pletely.It's not that I don't want to. It's that I've realizedthat I do not yet fully understand it myself :-]I mean, in a way that I can communicate toothers.I've been doing Maths in an analogous wayfor decades - just because it works - not re-alizing that there was this deeper stuff in thegeneral method I was using.I only became aware of that a couple of weeksago, and am still exploring it for myself.When I become better able to communicate itto others, I'll do that.[What I was actually doing in this discussionwas establishing a date of discovery.]> 9. This seems, to me, to be quite astounding,> because Incommensurability is nowhere to> be seen.> so, Continuously, without a hint of Incom-> mensurability.> are no longer a Problem - just construce> a Circle with the necessary radious, and, voil?,> there's your Conjugated-pair of Irrationals,> Numbers, but you can address them, Exactly, === Subject: Re: Yikes! I've Broken Incommensurability!In sci.math.num-analysis, ken:> In sci.math.num-analysis, ken> <4A7gc.13962$A_4.13209@newsread1.news.pas.earthlink.net>:> In sci.math.num-analysis, ken> :> [...]> In any event, if A=1 and B=1, how does C=sqrt(2)> become rational by moving the point on the semicircle> around, or mutating the radius?> [...]> It's in the Proof I posted.> Given a circle with radius sqr(2), because> sqr(2) = 4. Did you mean sqrt(2)? Admittedly> it depends on the language package; the only place> I've personally seen sqr() is in a variant of Pascal.sqr() is in BASIC.> A^2 + B^2 = C^2, always,> Given the constraints of the circle, yes.> [A[ca] + B[ca]] = sqr(2)> in which A[ca] and B[ca] are EquivalentIntegers =in the Number system that's> established in the Proof=, measurable in> the same Units of central angle, eachone unit, but positional.> It's what I'm 'excited ' about - a whole new> Number system, that's extremely-rich, and> in which all numbers are Integers is estab-> lished in the proof.> What establishes this Integer-ness is the> central angle [ca] that [for the purposes of this> discussion] is set, and used, unchanged, through-> out a calculation. All of the new Number sys-> tem's units are in terms of the central angle> that's employed.> It's a bit like Complex Numbers, A[ca] and> B[ca] are conjugated pairs, but it's all Inte-> gers[ca].> It's [another] whole new way of doing Maths.> Very flexible - because each central angle> establishes a different set of Integers.> And, obviously, above, they map Exactly to> traditional Integers.> If they map exactly then you might have a problem.> In the realm of standard integers:> Assume A/B = sqrt(2), and A and B are integers with no common factors> (i.e., no D > 1 such that A/D and B/D are integers).> Then A^2/B^2 = 2> and therefore A is even.> Write A = 2*C> then 4*C^2/B^2 = 2> or C^2/B^2 = 1/2> or B^2/C^2 = 2> and therefore B is even.> Therefore A and B have a common factor 2, contradicting> the original assertion that they don't.> Hence A/B = sqrt(2) is impossible for straight integers,> unless both A and B are equal to 0, which is fairly useless.> This [ca] Integer mapping may address this but if so it's> not clear to me exactly how.> Basically, it's why, as one scales simple sine-> and cosine-based plots smaller and smaller, no> discontinuities will ever show up, even though> there're standard difficult numbers all over the> place, in-there.> It's part of the bridge, between simple Geom-> etry and higher-level Maths, that I discussed> while discussing my Proof of Fermat's Last> Theorem.> What I envision is that, through bi-directional> transformations, these new Integers can be> Mapped into the 'sticky-parts' of all manner of> Problems, Exactly, thereby, enabling folks to> proceed beyond the 'sticky-parts'.> The problem above has a sticky-part -- namely, the> contradiction. Perhaps you can clarify with thisnew math how one resolves around this issue?> All I'm doing, at this 'time', is saying, Hey,> here's some new Integers. [I've been using> them for as long as I can remember, without> realizing that they were, in fact, Integers that> can, infact, transform to traditional Numbers.> In my Proof of Fermat's Last Theorem, for> instance, I Proved that, Given:> X^2 + Y^2 = Z^2> for all n > 2, either X[ca] or Y[ca] is al-> ways > 1 and [ca] is always > 0.]> I see no N and am not sure what X[ca] means. You'll> probably want to define precisely how one adds X[ca]> and Y[ca], subtracts, multiplies, and divides.> It's not a problem.Sure looks like one to me. :-)> Basically, all I'm doing is reducing myriad things> to the same central angle, and using the fact that> all manner of things can be reduced to the same> central angle, as a doorway through which one> can 'walk' to switch between other Maths.I can reduce all integers to the number 0 but thatwouldn't be a bijective mapping.> And it's quite a spectacular doorway> Of course, like any other technique, one must> expend some energy in constructing the door-> way Exactly, that's all [which, basically, consists> of casting Maths, on both sides of the doorway> as Fermat Numbers.> Then, voil.88! The Infinite-Exactness of 'the' Circle> allows one to do spectacular things in Maths.> [I did a version of the same stuff in my Proof of> Fermat's Last Theorem, only, in that Proof, us-> ing the Infinite-Exactness of Squares, instead of> Circles.> To get a hangle on the general methods that I'm> discussing, it's best to start with the FLT Proof,> which was discussed most-completely in> bionet.neuroscience [I'm an Amateur Neurosci-> entist].> Above, you're not casting the numbers into theFermat Integers, so, of course, it doesn't work.> [And, Forgive me, Please, I'm wanting to just> continue to reify the method.]> No arithmetical methods are any problem, be-> cause, all one has to do is carry-out the arith-> metic, getting an answer, which is, then, used> as the radius - voil.88! Fermat Integers, mapped> Exactly.> Last night, for the fun of it, I used:> ((SQR(2) ^ SQR(2)) ^ (EXP(1)) / pi)Need a ')'. :-)With that correction one gets, in standard math,1.5281947... Can't say if it's irrational ornot but it's definitely not an integer.> Voil.88! Fermat Integers mapped to any cen-> tral angle, which allows one to transform, via> the central angle.> What's neat about it is that it only requires two> more steps, +/- [ca]-variations, and everything> is Nailed-down Exactly, interms of relative lengths.> Which is what I was referring to when I said> that I'd Eliminated 'irrationals' - they just 'go away' -> lose their 'sting' - lose their quality of being able> to 'stymie' Solution.> Because everything in the method I've described> is mutually-Exact [which is why I've stated that> everything in the method is Integer].> I've been doing ZMaths this way for decades, but> it was only during the necessities inherent in Def-> ending my Proof of FLT that I explored more-> deeply into the 'guts' of the Method, and saw the> stuff that I'm discussing in this thread.> It's fully as big as The Calculus, and both, more-> powerful, and enormously-faster.> I'll discuss the Method further, in-person.You might want to put this on a webpage, and givea number of concrete examples with step-by-stepwalkthroughs on exactly how these problems arereduced -- and the correct answers.As it is, it appears to be a lot of handwaving, whichto me isn't all that useful.#191, ewill3@earthlink.netIt's still legal to go .sigless. === Subject: Re: Yikes! I've Broken Incommensurability!> In sci.math.num-analysis, ken> :> In sci.math.num-analysis, ken> <4A7gc.13962$A_4.13209@newsread1.news.pas.earthlink.net>:in> In sci.math.num-analysis, ken> :> [...]> In any event, if A=1 and B=1, how does C=sqrt(2)> become rational by moving the point on the semicircle> around, or mutating the radius?> [...]> It's in the Proof I posted.> Given a circle with radius sqr(2), because> sqr(2) = 4. Did you mean sqrt(2)? Admittedly> it depends on the language package; the only place> I've personally seen sqr() is in a variant of Pascal.sqr() is in BASIC.> A^2 + B^2 = C^2, always,> Given the constraints of the circle, yes.> [A[ca] + B[ca]] = sqr(2)> in which A[ca] and B[ca] are EquivalentIntegers =in the Number system that's> established in the Proof=, measurable in> the same Units of central angle, eachone unit, but positional.> It's what I'm 'excited ' about - a whole new> Number system, that's extremely-rich, and> in which all numbers are Integers is estab-> lished in the proof.> What establishes this Integer-ness is the> central angle [ca] that [for the purposes of this> discussion] is set, and used, unchanged, through-> out a calculation. All of the new Number sys-> tem's units are in terms of the central angle> that's employed.> It's a bit like Complex Numbers, A[ca] and> B[ca] are conjugated pairs, but it's all Inte-> gers[ca].> It's [another] whole new way of doing Maths.> Very flexible - because each central angle> establishes a different set of Integers.> And, obviously, above, they map Exactly to> traditional Integers.> If they map exactly then you might have a problem.> In the realm of standard integers:> Assume A/B = sqrt(2), and A and B are integers with no common factors> (i.e., no D > 1 such that A/D and B/D are integers).> Then A^2/B^2 = 2> and therefore A is even.> Write A = 2*C> then 4*C^2/B^2 = 2> or C^2/B^2 = 1/2> or B^2/C^2 = 2> and therefore B is even.> Therefore A and B have a common factor 2, contradicting> the original assertion that they don't.> Hence A/B = sqrt(2) is impossible for straight integers,> unless both A and B are equal to 0, which is fairly useless.> This [ca] Integer mapping may address this but if so it's> not clear to me exactly how.> Basically, it's why, as one scales simple sine-> and cosine-based plots smaller and smaller, no> discontinuities will ever show up, even though> there're standard difficult numbers all over the> place, in-there.> It's part of the bridge, between simple Geom-> etry and higher-level Maths, that I discussed> while discussing my Proof of Fermat's Last> Theorem.> What I envision is that, through bi-directional> transformations, these new Integers can be> Mapped into the 'sticky-parts' of all manner of> Problems, Exactly, thereby, enabling folks to> proceed beyond the 'sticky-parts'.> The problem above has a sticky-part -- namely, the> contradiction. Perhaps you can clarify with thisnew math how one resolves around this issue?> All I'm doing, at this 'time', is saying, Hey,> here's some new Integers. [I've been using> them for as long as I can remember, without> realizing that they were, in fact, Integers that> can, infact, transform to traditional Numbers.> In my Proof of Fermat's Last Theorem, for> instance, I Proved that, Given:> X^2 + Y^2 = Z^2> for all n > 2, either X[ca] or Y[ca] is al-> ways > 1 and [ca] is always > 0.]> I see no N and am not sure what X[ca] means. You'll> probably want to define precisely how one adds X[ca]> and Y[ca], subtracts, multiplies, and divides.> It's not a problem.> Sure looks like one to me. :-)> Basically, all I'm doing is reducing myriad things> to the same central angle, and using the fact that> all manner of things can be reduced to the same> central angle, as a doorway through which one> can 'walk' to switch between other Maths.> I can reduce all integers to the number 0 but that> wouldn't be a bijective mapping.> And it's quite a spectacular doorway> Of course, like any other technique, one must> expend some energy in constructing the door-> way Exactly, that's all [which, basically, consists> of casting Maths, on both sides of the doorway> as Fermat Numbers.> Then, voil.88! The Infinite-Exactness of 'the' Circle> allows one to do spectacular things in Maths.> [I did a version of the same stuff in my Proof of> Fermat's Last Theorem, only, in that Proof, us-> ing the Infinite-Exactness of Squares, instead of> Circles.> To get a hangle on the general methods that I'm> discussing, it's best to start with the FLT Proof,> which was discussed most-completely in> bionet.neuroscience [I'm an Amateur Neurosci-> entist].> Above, you're not casting the numbers into theFermat Integers, so, of course, it doesn't work.> [And, Forgive me, Please, I'm wanting to just> continue to reify the method.]> No arithmetical methods are any problem, be-> cause, all one has to do is carry-out the arith-> metic, getting an answer, which is, then, used> as the radius - voil.88! Fermat Integers, mapped> Exactly.> Last night, for the fun of it, I used:> ((SQR(2) ^ SQR(2)) ^ (EXP(1)) / pi)> Need a ')'. :-)> With that correction one gets, in standard math,> 1.5281947... Can't say if it's irrational or> not but it's definitely not an integer.> Voil.88! Fermat Integers mapped to any cen-> tral angle, which allows one to transform, via> the central angle.> What's neat about it is that it only requires two> more steps, +/- [ca]-variations, and everything> is Nailed-down Exactly, interms of relative lengths.> Which is what I was referring to when I said> that I'd Eliminated 'irrationals' - they just 'go away' -> lose their 'sting' - lose their quality of being able> to 'stymie' Solution.> Because everything in the method I've described> is mutually-Exact [which is why I've stated that> everything in the method is Integer].> I've been doing ZMaths this way for decades, but> it was only during the necessities inherent in Def-> ending my Proof of FLT that I explored more-> deeply into the 'guts' of the Method, and saw the> stuff that I'm discussing in this thread.> It's fully as big as The Calculus, and both, more-> powerful, and enormously-faster.> I'll discuss the Method further, in-person.> You might want to put this on a webpage, and give> a number of concrete examples with step-by-step> walkthroughs on exactly how these problems arereduced -- and the correct answers.> As it is, it appears to be a lot of handwaving, which> to me isn't all that useful.I Acknowledge that your Criticisms are Valid,because I've not, yet, stated things Completely. === Subject: Re: advanced numerical texts?> Can anyone suggest a text for topics covering implicit time-stepping ODE> methods? I have the text by Cheney and Klincaid and it doesn't cover it.> Specifically I'm interested in learning how to implement such a method on> -kIf you are looking for a application package to do this for you check out 1) xppaut : http://www.pitt.edu/~phase/2) dstool : http://www.cam.cornell.edu/guckenheimer/dstool.html3) content : ftp://ftp.cwi.nl/pub/CONTENTAll these are very stable codes used by reseachers in applied dynamicalsystems. They also happen to be free and open sourced with rpms and/ordebs available.As for a book my personal preferance is Butchers' Numerical methods forordinary differential equations. The books by E. Hairer and G. Wannerare ofcourse the bible. === Subject: Re: advanced numerical texts?There are many texts, and I don't have any one in particular to recommend.But, if you want to see such algorithms in action (step-by-step, doing theactual calculations on actual problems you enter -- including theroot-finding implicit part), you might take a look at a beta program I'mdeveloping on numerical ODEs, called NODE. It's not complete, but has a lotof stuff ready to use. You can get it from http://mason.gmu.edu/~bcrain/.Ben Crain> Can anyone suggest a text for topics covering implicit time-stepping ODE> methods? I have the text by Cheney and Klincaid and it doesn't cover it.> Specifically I'm interested in learning how to implement such a method on> -k === Subject: Re: advanced numerical texts?> Can anyone suggest a text for topics covering implicit time-stepping ODE> methods? I have the text by Cheney and Klincaid and it doesn't cover it.> Specifically I'm interested in learning how to implement such a method onMy personal favorite to learn the subject along with code is still thetextbook by C. W. Gear. I like compactness, and Gear presents keypoints in a few pages where other authors drone on forever , chapterafter chapter, and put you to sleep with irrelevant details dear totheir heart.The publication date (1971), however, means two problems: (1) out ofprint, used copies in net run $60-80, best bet is a library; (2) theillustration language is Fortran IV, which was then a logical choicebut now is a dead end. But the code is very well commented andtranslation to C, Java, Matlab, etc, should be straightforward. === Subject: need help on Normalizing dataI have a problem on working with infinite power series in Matlabi am using this seriesfn=summation(an*x^n)where x=50 and n=1,2,...........when x>182, the term x^n is going to infinity....This may be due tothe value is out of matlab limit (its a very big number)but I want to use more number of terms in my apllication (may bearound 350 terms)is it possible to overcome this problem by normalizing x??If yes, How can I do this in matlab?Please suggest me how to deal with this problemRavi === Subject: Re: need help on Normalizing data> I have a problem on working with infinite power series in Matlab> i am using this series> fn=summation(an*x^n)> where x=50 and n=1,2,...........> when x>182, the term x^n is going to infinity....This may be due to> the value is out of matlab limit (its a very big number)*wince*Isn't it just.What kind of an answer are you hoping to get out of this summation?Unless your {an} are pretty special, the result is almost certainlygoing to be meaninglessly large, and once rounding errors are takeninto account, any meaningful result will likely be completely swamped.If your x satisfied |x| < 1 it might work better. For x=50 I can'tsee it ever working. What is the application behind the numbers?Colin === Subject: bessel function with complex argument by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i3KGt2K15868; dear everone,Could you give me some information about bessel with complex argument? I need the C code! thank you! === Subject: Re: bessel function with complex argument > dear everone,Could you give me some information about bessel with >complex argument? I need the C code! > thank you! http://www.netlib.org search 'bessel'there are lots of entries (for example toms/644 fits your needs) however only cephes/besseltgz is in c. but this one has not the complex argument feature. so you smusttranslate toms/644 into c yourself which however should not be too hardmaybe f2c helps === Subject: C++ instead of DelphiThe main limitations of the CFD program I am developing are the number ofgrid points and the time it takes to run on a PC. Obviously the two arelinked as the run time goes up as number of points cubed roughly. The morepoints I can use the greater the accuracy. The main time in the program istaken up with evaluating influence coefficients of source and doublet panelsand then solving a set of simultaneous equations. Would there be anyadvantage in changing from programming in Delphi, which I am familiar with,to something like C++, which I would have to buy and learn? === Subject: Re: C++ instead of Delphi> The main limitations of the CFD program I am developing are the number of> grid points and the time it takes to run on a PC. Obviously the two are> linked as the run time goes up as number of points cubed roughly. The more> points I can use the greater the accuracy. The main time in the program is> taken up with evaluating influence coefficients of source and doublet panels> and then solving a set of simultaneous equations. Would there be any> advantage in changing from programming in Delphi, which I am familiar with,> to something like C++, which I would have to buy and learn?As you have hinted, the main problem seems to be algorithmic. The samealgorithms won't run much faster just because they are written inC++. Perhaps time is better spent finding better algorithms? Besides,I see no reason why Delphi should be much slower - ask an expert tohelp you with optimization.for g++), so you do not have to buy one.If you like Delphi, there is some chance that you will find C++awkward, ugly, and error prone, and may wonder why on earth people areso fond of it. I don't know an answer to that. === Subject: Re: C++ instead of Delphi> The main limitations of the CFD program I am developing are the number of> grid points and the time it takes to run on a PC. Obviously the two are> linked as the run time goes up as number of points cubed roughly. The more> points I can use the greater the accuracy. The main time in the program is> taken up with evaluating influence coefficients of source and doublet panels> and then solving a set of simultaneous equations. Would there be any> advantage in changing from programming in Delphi, which I am familiar with,> to something like C++, which I would have to buy and learn?> As you have hinted, the main problem seems to be algorithmic. The same> algorithms won't run much faster just because they are written in> C++. Perhaps time is better spent finding better algorithms? Besides,> I see no reason why Delphi should be much slower - ask an expert to> help you with optimization.> for g++), so you do not have to buy one.> If you like Delphi, there is some chance that you will find C++> awkward, ugly, and error prone, and may wonder why on earth people are> so fond of it. I don't know an answer to that.Neither do I. I like Delphi very much, for number-crunching as well as guibuilding. And I also find C++ awkward and ugly. About error prone I know not,since I refuse to use it. But I think the question has to do with the platform onwhich the program should run. If it is to be pc Windows, there is no betterchoice than Delphi. If this is too dogmatic a statement, I would appreciatehearing evidence to the contrary. === Subject: Re: C++ instead of Delphi>The main limitations of the CFD program I am developing are the numberof>grid points and the time it takes to run on a PC. Obviously the twoare>linked as the run time goes up as number of points cubed roughly. Themore>points I can use the greater the accuracy. The main time in theprogram is>taken up with evaluating influence coefficients of source and doubletpanels>and then solving a set of simultaneous equations. Would there be any>advantage in changing from programming in Delphi, which I am familiarwith,>to something like C++, which I would have to buy and learn?> As you have hinted, the main problem seems to be algorithmic. The same> algorithms won't run much faster just because they are written in> C++. Perhaps time is better spent finding better algorithms? Besides,> I see no reason why Delphi should be much slower - ask an expert to> help you with optimization.> for g++), so you do not have to buy one.> If you like Delphi, there is some chance that you will find C++> awkward, ugly, and error prone, and may wonder why on earth people are> so fond of it. I don't know an answer to that.> Neither do I. I like Delphi very much, for number-crunching as well asgui> building. And I also find C++ awkward and ugly. About error prone I knownot,> since I refuse to use it. But I think the question has to do with theplatform on> which the program should run. If it is to be pc Windows, there is nobetter> choice than Delphi. If this is too dogmatic a statement, I wouldappreciate> hearing evidence to the contrary.like it and other languages like C or its variants look much less easy tounderstand. The reason I asked is that in any discussion of which languageto use Delphi never gets mentioned and I was beginning to think there mightbe something wrong with it when everyone else seems to use C, C+, C++ orFortran. I used to program in Fortran on mainframes a long time ago butswitched to Turbo Pascal on a PC and followed its development into Delphi.I suppose the main decider for most programmers is that all the professionalnumerical analysis packages like netlib, lapack, or whatever they are, havebeen written in Fortran or a C variant. I have to get items likesimultaneous equation solvers out of Numerical Recipes. === Subject: Re: C++ instead of Delphi <4085CF3A.38902643@bellatlantic.net> like it and other languages like C or its variants look much less easy to> understand. The reason I asked is that in any discussion of which language> to use Delphi never gets mentioned and I was beginning to think there might> be something wrong with it when everyone else seems to use C, C+, C++ or> Fortran.What probably happens is that people who use other languages don'ttalk that much about it ;)> I suppose the main decider for most programmers is that all the professional> numerical analysis packages like netlib, lapack, or whatever they are, have> been written in Fortran or a C variant. I have to get items like> simultaneous equation solvers out of Numerical Recipes.Most languages have mechanisms for calling functions in otherlanguages. In particular, calling C and Fortran from Delphi should bepossible. Ask in comp.lang.delphi. === Subject: Re: Biharmonic problem by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i3KJ9HF01866;>In the biharmonic problem:>- boundary conditions on function and first derivative are Dirichlet>(essential or kinematic)>- boundary conditions on second and thid derivative are Von Neumann(natural>or static)>is that correct?>AntonioIn the biharmonic problem [b.p.] (contrary to the harmonic one)there is no division into Dirikhlet and Neumann (not John von Neumann, but Franz Neumanna German mathematician of 19th century)boundary conditions.The first one is sometimes called the classical b.p.In fact, the second one can be transformed into the first one.You may wish to look at my review paperSelected topics in the history of the b.p.Slava Meleshko === Subject: Re: Underdetermined system least squares>Hi folks,> I am doing some linear least square calculations. I have been using>the non-negative least squares algorithm written by Lawson and Hanson.>Unfortuantely I don't have a copy of their book, (its on order).>The problem involves seismic wave-form modeling. Depending on the>number of Green's functions I attempt to fit to the data, my problem>can either be overdetermined (if I use relatively few green's>functions) or underdetermined ( great gobs of green's functions).>Essentially I'm solving for the weights or amplitude (which can't be>negative) of the green's functions that result in the best fit to the>observed wave forms.>My question is do the solutions I get from the non-negative least>squares algorithm have any meaning when the system is underdetermined?>I'm not sure exactly what NNLS written by Lawson and Hanson do in that>case.>Truth be told, the results I get look reasonable. Course, that doesn't>mean they are not bull.>NNLS solves > ||Cx-c||=min subject to x>=0. > NNLS is an outdated method that takes O(n^4) work.> Just solving it as a bound constrained QP is much> faster once n is not very small.OK. Did not know that. Do you have a recommendation? I've heard ofLSSOL, but I haven't used it before. Do you think it would do thetrick?> result in a unique solution, since although this constraint isphysical, the green's functions are not the real Earth and hence> w/o the non-negative constraint there may well be x(i) < 0.> I was figuring that if the system was truly underdetermined, then the> result I get should be exact in the sense that the observed> waveforms can be precisely recovered using the determined> coefficients. I haven't seen that, therefore I conclude the system is> not truly underdetermined.> Typically, nonnegativity constraints regularize a little.> But if the system is significantly underdetermined and one> hasn't included regularization terms in your least squares> formulation, one can be almost sure to get artifacts, or at least> less resolution that could be achievable utilizing all the> qualitative information you have.OK. Well now I'll show off my ignorance. I don't know what aregularization term is. Is that some sort of smoothing? Well anyway,I will research it.By the way, I'd be interested if there is a way of minimizing thisobjective instead: [ 1. - (x.x)(b.b) / (b.b)(b.b) ] = where x is the solutionvector and b is the rhs. Is there a class of problems this falls under?> The image reconstruction paper in> http://www.mat.univie.ac.at/~neum/regul.html> uses nonnegativity constraints and discusses relevant issues,> though not for NNNLS but for a different method. I will have a look at it. Just to be sure I have the right one,this paper involved the conjugate gradient method?