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