solving nonlinear equation using newton method

2 views (last 30 days)
juwita14
juwita14 on 11 Dec 2021
Edited: Pavl M. on 20 Oct 2024
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

Answers (2)

John D'Errico
John D'Errico on 18 Oct 2024
Edited: John D'Errico on 18 Oct 2024
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.

Pavl M.
Pavl M. on 18 Oct 2024
Edited: Pavl M. on 20 Oct 2024

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
%+380990535261
%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

  2 Comments
Torsten
Torsten on 18 Oct 2024
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.
Pavl M. on 20 Oct 2024
Edited: Pavl M. on 20 Oct 2024
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.

Sign in to comment.

Categories

Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!