I am trying to solve a non linear simultaneous equation using NR method.The program does not converges.I know the exact solution and if the initial value is set to exact value then the program converges in one step.

 I have attached the m file with the initial values of the three unknowns equal to its exact values (i.e in the x vector) as a result the value converges in a single iteration.

 If i initialize value little different from the exact solution.Then the value does not converges.I want my program to converge with values different from exact initial value.The exact value Of X vector is

x=[0.033670275;0.251;1168]

If the initial X vector is changed from the above values then the program does not converges. Please help me in this regard.

Answer 1

[Are Mjaavatten]

First: Your script parameter_function_matrix.m is a bit confusing, with your exceedingly long expressions. It would be easier to read and debug if you defined some variables that you could use to build your function expressions. Example: For the first element of F you could write:

f = zeros(3,1);
a = 7.36;
b = 24.2;
c = 6.83;
d = 30.4;
f(1) = a-(b+x(2)*(c-a))/x(3)-(a-(d-a*x(2))/x(3))*exp((b+c*x(2)-d)/x(1)/50);

Convergence of the Newton-Raphson method may depend strongly on the initial values, but starting not too far from the known solution it converges in most cases. In the script below, I start at a random values normally distributed around the correct solution, using a standard deviation equal to the solution value. It converges in about half the cases, although sometimes to another solution. A narrower distribution will converge more often.

xsol = [0.033670275;0.251;1168];   % Sought solution
x0 = xsol.*(1+randn(3,1));  % Try starting a bit away from the solution

% Newton-Raphson:
x = x0;
for i = 1:50
    [f,J] = patro(x);
    dx = J\f(:);
    x = x - dx;
    if norm(f(:))< 1e-5
        break;
    end
end
if i >= 50
    disp('No convergence')
else
    fprintf('%12s %12s\n','x','g');
    for i = 1:3
        fprintf('%12g %12g\n',x(i),f(i));
    end
end

I have made some minimal edits to your (exhausting) function code, and attached it as Matlab function patro.m.

[sanat patro]

Sir the function matrix F,C and the J matrix explanation and details are given below

Below 3 equations are the actual set of equations

c1=((a1)-((a3+(a4*x(2))-(a1*x(2)))/(x(3)))-(((a1)-((a2-(a1*x(2)))/(x(3))))*exp((a3+(a4*x(2))-a2)/(a5*x(1)))));

c2=((((a3*B)*(exp(D)))+(a3/x(3)))/(1+((B*x(2))*exp(D))+(x(2)/x(3))));

c3=(((B*exp(E))+(1/x(3)))/(1+(B*x(2)*exp(E))+(x(2)/x(3))))*x(3);

where

D=((a3+(a4*x(2))-(a2))/(a5*x(1)));

B=(((a1*x(3))-a2+(a1*x(2)))/(a5*x(1)*x(3)));

E=(((a1*x(2))-a2)/(a5*x(1)));

F matrix contains the RHS of actual equation and C matrix contains the value of c1 ,c2 and c3.

c1=6.841128071170712 c2=6.828573650081349 c3=1.000433970354048

The J matrix is calculated as given below

%This is the way the jacobian matrix is calculated. %It is calculated by replacing x(1) with x1,x(2) with x2,and x(3) with x3. %Again after getting the J matrix x1,x2,x3 are replaced by x(1),x(2),and %x(3)

syms x1;

syms x2;

syms x3;

a1=7.36;

a2=30.4;

a3=24.2;

a4=6.83;

a5=50;

D=((a3+(a4*x2)-(a2))/(a5*x1));

B=(((a1*x3)-a2+(a1*x2))/(a5*x1*x3));

E=(((a1*x2)-a2)/(a5*x1));

f1=((a1)-((a3+(a4*x2)-(a1*x2))/(x3))-(((a1)-((a2-(a1*x2))/(x3)))*exp((a3+(a4*x2)-a2)/(a5*x1))));

f2=((((a3*B)*(exp(D)))+(a3/x3))/(1+((B*x2)*exp(D))+(x2/x3)));

f3=(((B*exp(E))+(1/x3))/(1+(B*x2*exp(E))+(x2/x3)))*x3;

j1=diff(f1,x1);

j2=diff(f1,x2);

j3=diff(f1,x3);

j4=diff(f2,x1);

j5=diff(f2,x2);

j6=diff(f2,x3);

j7=diff(f3,x1);

j8=diff(f3,x2);

j9=diff(f3,x3);

J=[j1 j2 j3 ;j4 j5 j6;j7 j8 j9];

And then The main program can be executed

I have again attached the main program

The exact value of x(1),x(2) and x(3) are

x(1)=0.033670275;

x(2)=0.251;

and x(3)=1168;

In the main program this is the value chosen and therefore the it takes single iteration to converge .But if we choose any other value near it the value does not converges.

The x values provided as

x=[0.033670275;0.251;1168]

• you please change the value of x as given below with bold in the mainprogram and you will notice that the program does not converges.These values are very near to the actual values

x=[0.5;0.5;1000] %Take this as initial estimate

And please suggest me how to converge.