solving nonlinear equation using newton method

hello can you help me with this

Function

f(x1,y1)4-x^2-y^2=0

f(x2,y2)1-e^x-y=0

Initial valuesx=1.00

y=-1.70

i am confused with the exponensial in the matlab code

%Function NewtonRaphson_nl() is given below.
fn = @(v) [(-v(1)^2)-v(2)^2+4 ; 1-v(4)^(x)-v(3);
    jacob_fn = @(v) [(-2*v(1)) 2*v(2); (-v(4)^(x)) -1; 
    error = 10^-5 ;
    
    
    v = [1 ;-1.7] ;
    no_itr = 20 ;
    
    [point,no_itr,error_out]=NewtonRaphson_nl(v,fn,jacob_fn,no_itr,error)
    NewtonRaphson_nl_print(v,fn,jacob_fn,no_itr,error);
    
    function [v1 , no_itr, norm1] = NewtonRaphson_nl(v,fn,jacob_fn,no_itr,error)
% nargin = no. of input arguments
if nargin <5 , no_itr = 20 ; end
if nargin <4 , error = 10^-5;no_itr = 20 ; end
if nargin <3 ,no_itr = 20;error = 10^-5; v = [1;1;1]; end

v1 = v;
fnv1 = feval(fn,v1);
i = 0;
while true
    jacob_fnv1 = feval(jacob_fn,v1);
    H = jacob_fnv1\fnv1;
    v1 = v1 - H;
    fnv1 = feval(fn,v1);
    i = i + 1 ;
    norm1 = norm(fnv1);
    if i > no_itr && norm1 < error, break , end
    %if norm(fnv1) < error , break , end
end
    end
    
    function [v1 , no_itr, norm1] = NewtonRaphson_nl_print(v,fn,jacob_fn,no_itr,error)
    v1 = v;
    fnv1 = feval(fn,v1);
    i = 0;
    fprintf('      Iteration|    x     |     y      |    z      | Error      | \n')
    while true
        norm1 = norm(fnv1);
        fprintf('%10d     |%10.4f| %10.4f | %10.4f| %10.4d |\n',i,v1(1),v1(2),v1(3),norm1)
        jacob_fnv1 = feval(jacob_fn,v1);
        H = jacob_fnv1\fnv1;
        v1 = v1 - H;
        fnv1 = feval(fn,v1);
        i = i + 1 ;
        norm1 = norm(fnv1);
        if i > no_itr && norm1 < error, break , end
        %if norm(fnv1) < error , break , end
        
    end
    end

Answer 1

[John D'Errico]

You just misunderstand how to use the exponential function.

In MATLAB, e^x is written as exp(x). That makes your fn as:

fn = @(v) [-v(1)^2-v(2)^2+4 ; 1-exp(v(1))-v(2)];

The jacobian is a MATRIX. A 2x2 matrix. The first row will be the partial derivatives of your first function, with respect to each variable. So the first column will be the derivatives with respect to v(1), the second column the derivatives with respect to v(2).

jacob_fn = @(v) [-2*v(1),-2*v(2);-exp(v(1)),-1];

So if you now pass in a vector v, TRY THEM OUT!

fn([1,2])

ans = 2×1

-1.0000 -3.7183

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

jacob_fn([1,2])

ans = 2×2

-2.0000 -4.0000 -2.7183 -1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Next, there is NO need to use feval. fn and jacob_fn are functions. Just pass in a vector of length 2, exactly as I did.

Finally, don't use variables named error. As then the FUNCTION named error will one day potentially fail. Choose your variable names wisely, as to not conflict with functions of that name.

Those were the obvious errors I saw, though I may have missed something else. It should get you close now.

The equations are quite simple, so we can even perform a symbolic solve. This way you will know if the Newton method has converged.

syms x y
[xsol,ysol] = vpasolve(fn([x,y]) == 0,[x,y],[1 1])
xsol =
-1.8162640688251505742443123715859
ysol =
0.8373677998912477276581914454592

Could there be other solutions? Of course. A nice trick is to use fimplicit. If the two curves cross, this is where a solution lives.

f12 = fn([x,y])
fimplicit(f12(1),'r')
hold on
fimplicit(f12(2),'b')
xlabel X
ylabel Y
hold off
grid on
axis equal

Ah, so there are indeed two solutions. And depending on your starting values chosen, you will find one or the other.

Answer 2

[Pavl M.]

I revamped the software code scientific computing TCE programme, here is good answer I provided:

clc
clear all
close all
function [v1 , no_itr, norm1] = NewtonRaphson_nl(v,fn,jacob_fn,no_itr,err)
% nargin = no. of input arguments
if nargin <5 , no_itr = 20 ; end
if nargin <4 , error = 10^-5;no_itr = 20 ; end
if nargin <3 ,no_itr = 20;error = 10^-5; v = [1;1;1]; end
v1 = v;
fnv1 = fn(v1);
i = 0;
while true
jacob_fnv1 = jacob_fn(v1);
H = jacob_fnv1\(fnv1+eps);
v1 = v1 - H;
fnv1 = fn(v1);
i = i + 1;
norm1 = norm(fnv1);
if i > no_itr || norm1 < err
break
end
%if norm(fnv1) < error , break , end
end
end

function [v1 , no_itr, norm1] = NewtonRaphson_nl_print(v,fn,jacob_fn,no_itr,err)
v1 = v;
fnv1 = fn(v1);
i = 0;
fprintf(' Iteration| x | y | Error | \n')
while true
norm1 = norm(fnv1);
fprintf('%10d |%10.4f| %10.4f | %10.4d | \n',i,v1(1),v1(2),norm1)
jacob_fnv1 = jacob_fn(v1);
H = jacob_fnv1\(fnv1+eps);
v1 = v1 - H;
fnv1 = fn(v1);
i = i + 1;
norm1 = norm(fnv1);
if i > no_itr || norm1 < err
break
end
%if norm(fnv1) < error , break , end
end
end
%Function NewtonRaphson_nl() is given below.
%For your provided input:
% f(x1,y1)4-x^2-y^2=0
% f(x2,y2)1-e^x-y=0
%I corrected for you:
fn = @(v)[(v(1)^2+v(2)^2)/4;v(2) + exp(v(1))];
jacob_fn = @(v)[v(1) v(2);exp(v(1)) 1];
err = 10^-7;
v = [1 ;-1.7];
no_itr = 25;
[point,no_itr,error_out]=NewtonRaphson_nl(v,fn,jacob_fn,no_itr,err)
[v1,no_itr,norm1] = NewtonRaphson_nl_print(v,fn,jacob_fn,no_itr,err)
%x = -0.1:0.01:0.1;
%y = x;

syms x y
[xsol,ysol] = vpasolve(fn([x,y]) == 0,[x,y],[1 1])

f12 = fn([x,y])
fimplicit(f12(1),'r')
hold on
fimplicit(f12(2),'b')
xlabel X
ylabel Y
hold off
grid on
axis equal
title('1st function')

%For 2nd function example:

fn = @(v) [-v(1)^2-v(2)^2+4 ; 1-exp(v(1))-v(2)];
jacob_fn = @(v) [-2*v(1),-2*v(2);-exp(v(1)),-1];

[point,no_itr,error_out]=NewtonRaphson_nl(v,fn,jacob_fn,no_itr,err)
[v1,no_itr,norm1] = NewtonRaphson_nl_print(v,fn,jacob_fn,no_itr,err)

syms x y
[xsol,ysol] = vpasolve(fn([x,y]) == 0,[x,y],[1 1])

f22 = fn([x,y])
fimplicit(f22(1),'r')
hold on
fimplicit(f22(2),'b')
xlabel X
ylabel Y
hold off
grid on
axis equal
title('2nd function')

fn = @(v)[2*v(1)^5+4;v(2) + exp(v(2))*exp(v(1))];

jacob_fn = @(v) [10*v(1)^4,0;exp(v(1))*exp(v(2)),exp(v(1))*exp(v(2))+1];

[point,no_itr,error_out]=NewtonRaphson_nl(v,fn,jacob_fn,no_itr,err)

[v1,no_itr,norm1] = NewtonRaphson_nl_print(v,fn,jacob_fn,no_itr,err)

syms x y
[xsol,ysol] = vpasolve(fn([x,y]) == 0,[x,y],[1 1])

f22 = fn([x,y])
fimplicit(f22(1),'r')
hold on
fimplicit(f22(2),'b')
xlabel X
ylabel Y
hold off
grid on
axis equal
title('3rd function')

%Constructed from needing help code by
%https://independent.academia.edu/PMazniker
%[phone]
%https://diag.net/u/u6r3ondjie0w0l8138bafm095b
%https://github.com/goodengineer
%https://orcid.org/0000-0001-8184-8166
%https://willwork781147312.wordpress.com/portfolio/cp/
%https://www.youtube.com/channel/UCC__7jMOAHak0MVkUFtmO-w
%https://nanohub.org/members/130066
%https://pangian.com/user/hiretoserve/
%https://substack.com/profile/191772642-paul-m

[Torsten]

Why do you work with a function that does not have a root ?

fn = @(v)[(v(1)^2+v(2)^2)/4;v(2) + exp(v(1))];

[Pavl M.]

I completed as per original function specified in the question post:

Function

f(x1,y1)4-x^2-y^2=0

f(x2,y2)1-e^x-y=0

Why does it matter?

• I have added recently solution with clear root.

You may substitute any other function into it:

For example:

fn = @(v)[2*v(1)^5+4;v(2) + exp(v(2))*exp(v(1))];

My solution will still work.