[Jprogramming] Non-scalar Newton-Raphson

[Don Kelly]

You are right on target - you need a matrix of derivatives  which is 
called a Jacobian matrix

http://mathworld.wolfram.com/Jacobian.html

In the case of an example:  that  of a power system flow of real and 
reactive power as related to voltage phase and magnitude The Newton 
Raphson approach involves the Jacobian matrix.
This approach has been used since the late 1950's  for systems with up 
to000 to 10000 variables (something that one numerical analysis text 
dismissed, erroneously,  as not being  a realistic problem ) and 
provides fast convergence (partially due to the sin/cos differential cycle.
I am not a mathematician but an engineer using this as a tool. I have 
written this in Turbo Basic as well as in APL (STSC version11-no  
complex number capability). The APL version took about 11 lines of 
actual solution with added lines of checking limits and dressing up 
output. The other version involves about 3 pages of code.  I have 
partially developed a J version but overall I may not have the skill to 
make it more compact than the APL coding.

https://en.wikipedia.org/wiki/Power-flow_study

Don Kelly


On 2/13/2016 8:00 AM, Roger Hui wrote:
> If the atoms of the function are independent, then it's easiest to do f"0
> Newton (or f Newton"0) and be done with it.  If they are not independent,
> that is f is what's called a vector function in conventional mathematical
> notation, then you need to use %. (matrix divide) instead of % (plain old
> divide), among other things.
>
> My knowledge of vector calculus is minimal.  Perhaps the mathematicians in
> this Forum can help you if you really do have a vector function.
>
>
>
>
> On Sat, Feb 13, 2016 at 7:50 AM, Louis de Forcrand <olius at bluewin.ch> wrote:
>
>> I’ve been trying to write a conjunction that will find the zeros to a
>> function using the
>> Newton-Raphson method. The simplest way to do this is probably:
>>
>> English:
>> x_n+1 = x_n - f ( x_n ) / f ‘ ( x_n )
>>
>> J:
>> eznewt=: 2 : ‘ ( - u % u D.1 ) ^: n ‘
>>
>> This works fine for scalar -> scalar functions, but won’t work if the rank
>> of the
>> result is AND the rank of the arguments are above 1.
>> Probably the most evident situation where this would be a problem is if one
>> were searching for a minimum instead of a zero, in which case the algorithm
>> would be applied to the derivative of the function: +/@:*: D.1 newt 30 ] 3
>> 4
>> As I understand it, this would give a result shape at each iteration (with
>> a function u) of:
>>
>> ( $y ) , $ u eznewt 1 y
>>
>> where A is the original argument, or the result of the previous iteration.
>> What we would
>> want is a result with shape $y .
>>
>> First, let’s get some rules clear:
>> - the syntax should be: (verb) newt (# of iterations) (argument)
>> - the argument can be of any rank
>> - the shape of the argument matches the shape of the result at any
>> iteration
>> - this implies that $ u y matches i. 0 or $ y
>>
>> The hard part is getting u % u D.1 to have the shape of the argument.
>> If y is a vector, and u y is a scalar, then u D.1 y will be a vector, and
>> u D.1 D.1 y will be a $y by $y matrix. All that is needed then is to take
>> its diagonal with (<0 1)&|: .
>>
>> But what if y is a matrix? Since $ u D.1 y -: ( $ y ) , $ u y , I though
>> maybe running
>> y's axis together with |: might work
>>
>> newt=: 2 : '(- u diag_and_divide u D.1)^:n’
>> diag_and_divide=: [ % ] |:~ -@#@$@[ <@}. i.@#@$@]
>>
>>   but something about it doesn’t.
>> My head is $pinning now, and I figured I’d send this and then take a break.
>>
>> Thanks in advance!
>> Louis
>> ----------------------------------------------------------------------
>> For information about J forums see http://www.jsoftware.com/forums.htm
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm