solving nonlinear equation using newton method

5 次查看（过去 30 天）
juwita14
juwita14 2021-12-11
编辑： Pavl M. 2024-10-20
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

请先登录，再回答此问题。

回答（2 个）

John D'Errico
John D'Errico 2024-10-18
编辑：John D'Errico 2024-10-18
Ran in:
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. 2024-10-18
编辑：Pavl M. 2024-10-20

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 个评论
Torsten
Torsten 2024-10-18
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. 2024-10-20
编辑：Pavl M. 2024-10-20
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by