2

I tried non-linear polynomial functions and this code works well. But for this one I tried several methods to solve the linear equation df0*X=f0 using backslash or bicg or lsqr, also tried several initial values but the result never converge. 

% Define the given function
syms x1 x2 x3

x=[x1,x2,x3];

f(x)=[3*x1-cos(x2*x3)-1/2;x1^2+81*(x2+0.1)^2-sin(x3)+1.06;...
    exp(-x1*x2)+20*x3+1/3*(10*pi-3)];

% Define the stopping criteria based on Nither or relative errors

tol=10^-5; 
Niter=100;

df=jacobian(f,x);

x0=[0.1;0.1;-0.1];

% Setting starting values

error=1; 
i=0; 

% Start the Newton-Raphson Iteration

while(abs(error)>tol)

f0=eval(f(x0(1),x0(2),x0(3)));

df0=eval(df(x0(1),x0(2),x0(3))); 

xnew=x0-df0\f0; % also tried lsqr(df0,f0),bicg(df0,f0)

error=norm(xnew-x0);

x0=xnew;

i=i+1

if i>=Niter

    fprintf('Iteration times spill over Niter\n');

    return;

end

end
Jarvis
  • 35
  • 5
  • Plot the function and see about a better choice for a starting guess. – duffymo Feb 21 '17 at 20:24
  • the function is a vector comprised of 3 equations, how to plot ? – Jarvis Feb 21 '17 at 20:45
  • Why are you using `eval` there?!? It's completely unnecessary and will destroy your computation speed, as well as drawn you into the abysmal pit that is `eval` and its associated bad habits. – Adriaan Feb 21 '17 at 20:46
  • If I do not use eval it will be super slow and even stuck – Jarvis Feb 22 '17 at 01:00

1 Answers1

1

You'll need anonymous functions here to better do the job (we mentioned it in passing today!).

First, let's get the function definition down. Anonymous functions are nice ways for you to call things in a manner similar to mathematical functions. For example,

f = @(x) x^2;

is a squaring function. To evaluate it, just write like you would on paper f(2) say. Since you have a multivariate function, you'll need to vectorize the definition as follows:

f(x) = @(x) [3*x(1) - cos(x(2) * x(3)) - 1/2; ...

For your Jacobian, you'll need to use another anonymous function (maybe call it grad_f) and compute it on paper, then code it in. The function jacobian uses finite differences and so the errors may pile up with the Jacobian is not stable in some regions.

The key is to just be careful and use some good coding practices. See this document for more info on anonymous functions and other good MATLAB practices.

Sean Roberson
  • 142
  • 7
  • Thank you so much! – Jarvis Feb 23 '17 at 01:00
  • I still have a question, the gradient of f(x) I computed is the same with the result of jacobian function. Also, I only used jacobian in the first beginning not in the loop, will that count? – Jarvis Feb 23 '17 at 03:46