Thread Subject: Newton Raphson for nxn system using symbolic toolbox

Hello

I have already implemented the Newton Raphson method for a
system of n equations. However in my implementation I have
to give the function and its Jacobian (written as a
function) as a input variable.

Now I want to use the symbolic toolbox in order not to
calculate the Jacobian myself. My basic problem is later to
evaluate that expression, since I don't know before hand the
size of the system. While this is not a problem for an
"ordinary" function which uses vectors it seems to be a
problem for symbolic expression.


Here is an example trying to use vectors
 clear all
 syms x y
 e=[x y]
 f1=sin(e(1)-e(2))
 e(1)=1
 e(2)=-1
 eval(f1)

But this does not work! It gives sin(x+y)

While of course
 x=1
 y=-1
 eval(f1)
works and gives
f1= 0.9093
but of course for a nxn system I don't know before hand the
size of the system and how many variables I have to set to
numerical values before evaluating the expression f1.

Can anybody help me?

thanks

Uwe Brauer

On 11/1/2011 6:20 AM, Uwe Brauer wrote:
> Hello
>
> I have already implemented the Newton Raphson method for a
> system of n equations. However in my implementation I have
> to give the function and its Jacobian (written as a
> function) as a input variable.
>
> Now I want to use the symbolic toolbox in order not to
> calculate the Jacobian myself. My basic problem is later to
> evaluate that expression, since I don't know before hand the
> size of the system. While this is not a problem for an
> "ordinary" function which uses vectors it seems to be a
> problem for symbolic expression.
>
>
> Here is an example trying to use vectors
> clear all
> syms x y e=[x y]
> f1=sin(e(1)-e(2)) e(1)=1
> e(2)=-1
> eval(f1)
>
> But this does not work! It gives sin(x+y)
>
> While of course
> x=1
> y=-1
> eval(f1)
> works and gives f1= 0.9093
> but of course for a nxn system I don't know before hand the
> size of the system and how many variables I have to set to
> numerical values before evaluating the expression f1.
>
> Can anybody help me?
>
> thanks
>
> Uwe Brauer

Perhaps this Optimization Toolbox documentation example will help:
http://www.mathworks.com/help/toolbox/optim/ug/brn4nh7.html#brv_i_1

Alan Weiss
MATLAB mathematical toolbox documentation

On 01.11.11 11:20, Uwe Brauer wrote:

> Here is an example trying to use vectors
> clear all
> syms x y
> e=[x y]
> f1=sin(e(1)-e(2))
> e(1)=1
> e(2)=-1
> eval(f1)

Instead of assigning, try

subs(f1, e, [1 -1])

Or, as Alan pointed out, convert f1 to a MATLAB function. If you have
multiple evaluations of the same symbolic expression, that will usually
be faster.


Christopher

Christopher Creutzig <Christopher.Creutzig@mathworks.com> wrote in message <4EBA8B5A.5090609@mathworks.com>...
> On 01.11.11 11:20, Uwe Brauer wrote:
>
> > Here is an example trying to use vectors
> > clear all
> > syms x y
> > e=[x y]
> > f1=sin(e(1)-e(2))
> > e(1)=1
> > e(2)=-1
> > eval(f1)
>
> Instead of assigning, try
>
> subs(f1, e, [1 -1])
>
> Or, as Alan pointed out, convert f1 to a MATLAB function. If you have
> multiple evaluations of the same symbolic expression, that will usually
> be faster.
>
>
> Christopher

unfortunately my Matlab version does not provide this feature
 for symbolic expressions

But my main problem remains.

I can perfectly write matlab code for solving Newton Raphson for any nxn system, since Matlab uses dynamic memory,
and I have to use as a entry variable the system I want to solve and its jacobian.

However when i only want to have the system and use Matlabs symbolic tool to calculate the jacobian I do not know how to do that since I don't know the size of the jacobian in advance.

Uwe

On 11.11.11 16:17, Uwe Brauer wrote:

> I can perfectly write matlab code for solving Newton Raphson for any nxn system, since Matlab uses dynamic memory,
> and I have to use as a entry variable the system I want to solve and its jacobian.
>
> However when i only want to have the system and use Matlabs symbolic tool to calculate the jacobian I do not know how to do that since I don't know the size of the jacobian in advance.

I still don't understand the problem. Something like the following
sketch should work, no?

...
J = jacobian(expression, vars);
...
% within some loop:
  ...
  Fval = double(subs(expression, vars, values));
  Jac = double(subs(J, vars, values));
  ...



Christopher

Christopher Creutzig <Christopher.Creutzig@mathworks.com> wrote in message <4EC235EE.9040108@mathworks.com>...
> On 11.11.11 16:17, Uwe Brauer wrote:
>
> > I can perfectly write matlab code for solving Newton Raphson for any nxn system, since Matlab uses dynamic memory,
> > and I have to use as a entry variable the system I want to solve and its jacobian.
> >
> > However when i only want to have the system and use Matlabs symbolic tool to calculate the jacobian I do not know how to do that since I don't know the size of the jacobian in advance.
>
> I still don't understand the problem. Something like the following
> sketch should work, no?
>
> ...
> J = jacobian(expression, vars);
> ...
> % within some loop:
> ...
> Fval = double(subs(expression, vars, values));
> Jac = double(subs(J, vars, values));
> ...
>
>
>
> Christopher

right, after thinking about it I came up with a similar solution,
BTW why do you use double?

function xvec=newtoniteratesymb(f,e,x0,TOL,nmax)

n=x0;
fev=subs(f,e,n);
jac=jacobian(f,e);
jacev=subs(jac,e,n);
xvec=x0;
k=1;
diff=100;
 
while diff>TOL & k<nmax
  fev=subs(f,e,xvec(:,k));
  jacev=subs(jac,e,xvec(:,k));
  xvec(:,k+1)=xvec(:,k)- jacev\fev;
  diff=max(max(xvec(:,k+1)-xvec(:,k)));
  k=k+1;
end

thanks

Uwe

"Uwe Brauer" wrote in message <j9te3u$cqg$1@newscl01ah.mathworks.com>...
> Christopher Creutzig <Christopher.Creutzig@mathworks.com> wrote in message <4EC235EE.9040108@mathworks.com>...
> > On 11.11.11 16:17, Uwe Brauer wrote:
> >
> > > I can perfectly write matlab code for solving Newton Raphson for any nxn system, since Matlab uses dynamic memory,
> > > and I have to use as a entry variable the system I want to solve and its jacobian.
> > >
> > > However when i only want to have the system and use Matlabs symbolic tool to calculate the jacobian I do not know how to do that since I don't know the size of the jacobian in advance.
> >
> > I still don't understand the problem. Something like the following
> > sketch should work, no?
> >
> > ...
> > J = jacobian(expression, vars);
> > ...
> > % within some loop:
> > ...
> > Fval = double(subs(expression, vars, values));
> > Jac = double(subs(J, vars, values));
> > ...
> >
> >
> >
> > Christopher
>
> right, after thinking about it I came up with a similar solution,
> BTW why do you use double?
>
> function xvec=newtoniteratesymb(f,e,x0,TOL,nmax)
>
> n=x0;
> fev=subs(f,e,n);
> jac=jacobian(f,e);
> jacev=subs(jac,e,n);
> xvec=x0;
> k=1;
> diff=100;
>
> while diff>TOL & k<nmax
> fev=subs(f,e,xvec(:,k));
> jacev=subs(jac,e,xvec(:,k));
> xvec(:,k+1)=xvec(:,k)- jacev\fev;
> diff=max(max(xvec(:,k+1)-xvec(:,k)));
> k=k+1;
> end
>
> thanks
>
> Uwe

Hello, I head a similar problem some days ago and i have to say this method works perfectly for small systems but when you go after bigger problems you might find some problems while using the function jacobian(f) and computing it's value every iteration. This is very demending on memory and CPU.

On 15.11.11 11:15, Uwe Brauer wrote:
> BTW why do you use double?

To make sure the result is not a sym object. Often, subs will do that
when substituting double variables, but only if evaluation within MATLAB
works – in which case the double command does no harm either. If the
substitution needs to go through the symbolic toolbox engine, the result
will be a sym object:

>> syms x
>> subs(erf(i*x), x, 1)

ans =

erf(i)

>> double(subs(erf(i*x), x, 1))

ans =

                  0 + 1.650425758797543i


Christopher