Re: Re: Re: Nestwhile IHB

*To*: mathgroup at smc.vnet.net*Subject*: [mg57438] Re: [mg57392] Re: [mg57229] Re: [mg57201] Nestwhile IHB*From*: Pratik Desai <pdesai1 at umbc.edu>*Date*: Sat, 28 May 2005 05:39:00 -0400 (EDT)*References*: <EC9C79ED705B7D498CE6892836AA3CAA5A2340@EAHQ-CSMB9.rws.ad.ea.com> <200505260831.EAA18752@smc.vnet.net> <opsrd6nux9iz9bcq@monster.ma.dl.cox.net> <4296005F.60904@umbc.edu> <opsrehzyk9iz9bcq@monster.ma.dl.cox.net> <opsrey42c0iz9bcq@monster.ma.dl.cox.net>*Sender*: owner-wri-mathgroup at wolfram.com

DrBob wrote: > On my point #10, there's something wrong. It seems that > rep1[dw]/D[rep1[dw]] == 1, but that can't be true unless rep1 is the > exponential function E^dw. Maybe this is a Mathematica bug! (In my > version 5.1.1, anyway.) > > I think you're trying to use the Newton method, but that would be > something like > > Clear[dw,func] > func[dw_] = dw - rep1[dw]/rep1'[dw] > dw = FixedPoint[func, 0.05] > > Unfortunately, that leads to a crash at my machine. > > This does a little better: > > dw =. ; > startingValue=10.^(-18) > func[dw_] := (Print[dw]; > Re[dw - rep1[dw]/rep1'[dw]]) > Block[{$MinPrecision = 20, $MaxPrecision = 20}, > dw = FixedPoint[func, startingValue]] > > but it gets close to zero for some starting values without actually > converging, and may not converge at all for other starting values. I > doubt that 0 is the solution you're looking for. > > Mathematica has solvers that will undoubtedly do better than your > overt Newton's method, but I'm not certain what equation you're trying > to solve. > > Bobby > > On Thu, 26 May 2005 16:01:00 -0500, DrBob <drbob at bigfoot.com> wrote: > >> 1) I don't have the VisualDSolve package, so my mileage may vary. I >> don't think it mattered, though, so maybe you're not actually using >> the package? >> >> 2) You're still sorting twice in defining c, I notice. The second >> sort undoes the first, so eliminate one of them (the first, if you >> don't want the result to change). >> >> 3) In defining Ao, it looks as if you intend j to be a variable, but >> that definition doesn't actually work, if j didn't already have a >> value, or if you changed it later. You could change "=" to ":=" to >> fix that problem. This would be a clearer definition, too: >> >> Ao[j_, dw_] := PadRight[{dw}, 2j + 2] >> >> 4) Definitions like this make no sense: >> >> cp[j_] = cp[j] = D[c[j], t]; >> >> That, again, seems to treat j as a variable but doesn't work if it >> actually VARIES. >> >> What you may mean is this: >> >> cp[j_] := cp[j] = D[c[j], t] >> >> 5) Everywhere you've used N, it doesn't do anything. >> >> 6) To speed things up, eliminate ":=" where you can. Like this, for >> instance: >> >> r[dw_] = Integrate[c[j]*(wo^2*uopp[j] + >> uo[j] - e*wo*(1 - uo[j]^2)*uop[j]), >> {t, 0, 2*Pi}] >> >> and >> >> q[dw_] = Integrate[ >> c[j]*(e*(1 - uo[j]^2)*uop[j] - 2*wo*uopp[j]), >> {t, 0, 2*Pi}] >> >> and >> >> k[dw_] = Integrate[Outer[Times, c[j], >> wo^2*cpp[j] - e*wo*(1 - uo[j]^2)*cp[j] + >> (1 + 2*e*wo*uo[j]*uop[j])*c[j]], >> {t, 0, 2*Pi}] >> >> That does the integration ONCE rather than again every time you >> change j or dw... and it's a SLOW integration, so this really >> matters. This requires j to never change, but you obviously aren't >> changing it anyway... except by starting all over. >> >> 7) With your original code, in this expression: >> >> sol1 = Solve[expr1[j, Î?w, dw], coeff[j]] // Flatten // Sort >> >> Solve returns one solution, with ten Rules in a sublist. Flattening >> gets rid of the sublist, but so would >> >> sol1 = First@Solve[expr1[j, Î?w, dw], coeff[j]] >> >> Sorting makes no difference, since anything/.sol1 comes out the same, >> no matter what order the variables are in. Try to get out of the >> habit of using //Flatten//Sort for EVERYTHING, whether it needs it or >> not. >> >> 8) With number (6) in mind, it's probably best to eliminate j as a >> variable entirely. It's also best to state that dw IS a variable, >> wherever it appears. See attachment. >> >> 9) re[dw] and deltaAn[dw] involve powers of dw from one to 30, and >> this could be very unstable, numerically. rep1 is the square root of >> a 60-degree polynomial. So beware. >> >> 10) func1 doesn't depend on dw, because things cancel out. It's a >> constant (see attachment), so your FixedPoint statement doesn't do >> anything but set dw equal to that constant. >> >> 11) Once you've set dw equal to something, this statement: >> >> u[dw_, t_] = Total[c.Ao[dw] + c.Î?An[dw]] >> >> pretends to treat dw as a variable, but actually doesn't. Evaluate it >> without a semicolon at the end, and you'll see there's no dependence >> on dw. Use ":=" if you'll want dw to vary after this, or else >> eliminate dw_ on the left. (See attachment.) >> >> 12) I added a new-line to the PlotLabel. >> >> Bobby >> >> On Thu, 26 May 2005 12:59:11 -0400, Pratik Desai <pdesai1 at umbc.edu> >> wrote: >> >>> DrBob wrote: >>> >>>> The HTML version doesn't do us any good. >>>> >>>> Bobby >>>> >>>> On Thu, 26 May 2005 04:31:47 -0400 (EDT), Pratik Desai >>>> <pdesai1 at umbc.edu> wrote: >>>> >>>>> Barthelet, Luc wrote: >>>>> >>>>>> How is it going? Do you still need help? >>>>>> Dr Bob pointed out a few things, so were you able to fix them? >>>>>> If so you should upload the new notebook. >>>>>> >>>>>> You might also want to add as much english text as possible to >>>>>> explain >>>>>> what you are trying to do. Maybe even line per line. I still do not >>>>>> understand what you are trying to do and what results you would >>>>>> like to >>>>>> see. >>>>>> >>>>>> How are you going to know you have the proper result? >>>>>> Maybe we should study a bit that function d[e,dw]? Do you know >>>>>> what it >>>>>> is supposed to look like? >>>>>> >>>>>> Some of the stuff. >>>>>> >>>>>> C[j_] has a t variable in there is that normal? >>>>>> Same for uo[j_] >>>>>> The delta is using delta and deltab which are undefined. >>>>>> The integrals are using p as a variable that is not defined. >>>>>> Plus the other issues that Dr Bob mentioned of course. >>>>>> >>>>>> If you stick to this we will help you get to a successful result. >>>>>> >>>>>> Cheers, >>>>>> Luc >>>>>> >>>>>> >>>>>> >>>>>> >>>>> Dr Bob and Luc >>>>> Well here we go again >>>>> You guys were right Nestwhile was not needed, although the idea of >>>>> Nestuntil is not an entirely bad one. Anyhow, I was partially >>>>> successful >>>>> in finishing the project, if you are interested please click on the >>>>> image or nb file at wiki >>>>> >>>>> http://www.mathematica-users.org/webMathematica/wiki/wiki.jsp?pageName=Notebook:ihb2.nb >>>>> >>>>> >>>>> >>>>> I am getting closer to what I want but my computer crashes every >>>>> time I >>>>> try to increase the number of terms. I can only do up to j=4 i.e 9 >>>>> terms . >>>>> >>>>> Any takers :-) >>>>> >>>>> Thanks again for all of your help >>>>> >>>>> Until next time >>>>> >>>>> Pratik >>>>> >>>> >>>> >>>> >>> Dr Bob >>> >>> Here is the link for the nb file >>> >>> http://www.mathematica-users.org/mediawiki/images/c/cc/ihb2.nb >>> >>> Best Regards >>> >>> Pratik >>> >>> >> >> >> > > > Dr Bob Thanks again for all the help..FYI, IHB stands for Incremental Harmonic Balance, it is a part numerical and part analytical method to solve nonlinear diffferential equation like Duffing, Van der Pol etc.. Essentially you start out with your ode then you apply a pertubation on the variables like (Newton-Raphson Method) u= uo+delu (displacement) w=wo+del Then apply the gallerkin procedure to get a set of algebraic equation. Then find a numeric solutions for the dels by assuming certain initial conditions, The method was popularized in the 80's but as you can probably tell NDSolve is much quicker. > but it gets close to zero for some starting values without actually > converging, and may not converge at all for other starting values. I > doubt that 0 is the solution you're looking for. You are right, the idea is that the tolerance for Norm[re[dw]]does not change, that is why I used fixed point. I guess what we are looking for is the local minimum not the global one. As a result your initial value of 10^-18 give zero, the trivial solution If you are interested I can suggest some papers for further information Thanks again Pratik Desai -- Pratik Desai Graduate Student UMBC Department of Mechanical Engineering Phone: 410 455 8134

**References**:**Re: Re: Nestwhile IHB***From:*Pratik Desai <pdesai1@umbc.edu>