ODE with Newton method

27 views (last 30 days)
Roberto Facundo Luna Mallea
Commented: Torsten on 15 Sep 2022
I am trying to solve the ODE y'=x*cos(x^2)*y^2 using newton's method but my code keeps running forever.
Please, I really appreciate some help with this.
clear all
close all
n=1001
y0=1;
xmin=0; xmax=10;
x=linspace(xmin,xmax,n); %grid
step=(xmax-xmin)/n; %Step Size
yy=zeros(1,n); %Array for results for EEE
yy(1) = y0; % Initial value
%Define the Function and derivative
f=@(x,y) x*cos(x^2)*y^2;
df=@(x,y) 2*x*cos(x^2)*y ;
F=@(x,y,yn) y-yn-step*f(x,y); %Function to set to zero
J=@(x,y) 1 - step*df(x,y); %Jacobian
tol=1.e-5; % Tolerance
for i=xmin+1:n-1
yy(i+1)=yy(i); %initial guess for Newton's method
res=-J(yy(i+1))\F(yy(i+1),yy(i));
while (norm(res,inf)>1.e-10)
yy(i+1)=yy(i+1) + res;
res=-J(yy(1,1))\F(yy(i+1),yy(i));
end
yy(i+1)= yy(i+1) + res;
end
plot(x,yy,'k--')
xlabel('t')
ylabel('y')
Thanks
  1 Comment
Torsten
Torsten on 15 Sep 2022
For comparison:
fun = @(x,y)x*cos(x^2)*y^2;
n = 1001;
xmin=0; xmax=10;
x=linspace(xmin,xmax,n); %grid
y0 = 1;
[X,Y] = ode45(fun,x,y0);
plot(X,Y)

Sign in to comment.

Answers (1)

VBBV
VBBV on 15 Sep 2022
Edited: VBBV on 15 Sep 2022
clear all
close all
n=1001
n = 1001
y0=1;
xmin=0; xmax=10;
x=linspace(xmin,xmax,n); %grid
step=(xmax-xmin)/n; %Step Size
yy=zeros(1,n); %Array for results for EEE
yy(1) = y0; % Initial value
%Define the Function and derivative
f=@(x,y) x*cos(x^2)*y^2;
df=@(x,y) 2*x*cos(x^2)*y ;
F=@(x,y,yn) y-yn-step*f(x,y) %Function to set to zero
F = function_handle with value:
@(x,y,yn)y-yn-step*f(x,y)
J=@(x,y) 1 - step*df(x,y) %Jacobian
J = function_handle with value:
@(x,y)1-step*df(x,y)
tol=1.e-5; % Tolerance
for i=xmin+1:n-1
yy(i+1)=yy(i); %initial guess for Newton's method
res=-J(yy(i+1),x(i))\F(x(i),yy(i+1),yy(i)); % why missing input arguments ???
while (norm(res,inf)>1.e-10)
yy(i+1)=yy(i+1) + res;
res=-J(yy(1,1),x(i))\F(x(i),yy(i+1),yy(i));
end
yy(i+1)= yy(i+1) + res;
end
plot(x,yy,'k--')
xlabel('t')
ylabel('y')

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!