ODE with Newton method

13 次查看（过去 30 天）
Roberto Facundo Luna Mallea
评论： Torsten 2022-9-15
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 个评论
Torsten
Torsten 2022-9-15
Ran in:
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)

请先登录，再进行评论。

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

回答（1 个）

VBBV
VBBV 2022-9-15
编辑：VBBV 2022-9-15
Ran in:
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')

类别

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

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by