Using Implicit Euler Method with Newton-Raphson method

2022 5 11
1 个回答
回答已采纳
更新时间：2022 11 20
21 次查看（30 天）

1 个投票

So I'm following this algorithm to write a code on implicit euler method
and here is my attempt
function y = imp_euler(f,f_y,t0,T,y0,h,tol,N)
t = t0:h:T;
n = length(t);
y = zeros(n,1);
y(1) = y0;
for k = 1:n-1
g = @(z) z - y(k) - h*f(t(k+1),z);
gp = @(z) 1 - h*f_y(t(k+1),z) ;
y(k+1) = newton(f,f_y,y(k),tol,N);
end
end
where
function sol=newton(f,fp,x0,tol,N)
i=0;
sol = zeros(N,2);
fc=abs(f(x0));
while (fc>tol)
xc = x0 - (f(x0)/fp(x0));
fc=abs(f(xc));
x0 = xc;
i=i+1;
sol(i,:) = [i; x0];
if (i>N)
fprintf('Method failed after %d iterations. \n',N);
break
end
end
sol = sol(any(sol,2),:);
end
Unfortunately, my code does not work for some reason. Could anybody guide me on how to code this? Comments are appreciated.

7 个评论
显示 5更早的评论 隐藏 5更早的评论

Torsten
Torsten 2022-5-11
You forgot to include the call to imp_euler.
Torsten
Torsten 2022-5-11
I mean: The code you provided doesn't work without the knowledge of f, fy, t0,T,y0,h,tol and N (I think n).
Torsten
Torsten 2022-5-11
How did you come the conclusion that the code does not work if you didn't test it ? Because for testing, inputs of the mentionned variables had to be provided.
Eugene Miller
Eugene Miller 2022-5-11
i tested it with this
f = @(t, y) -20*t*y^2;
f_y = @(t, y) -40*t*y;
t0 = 0;
T = 1;
y0 = 1;
h = 0.2;
tol = 1e-8;
N = 100;
Eugene Miller
Eugene Miller 2022-5-11
and the answer is
y =
0.2
(but it should be a vector)

请先登录，再进行评论。

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

 采纳的回答

Torsten
Torsten 2022-5-11
编辑：Torsten 2022-5-11

1 个投票

f = @(t,y) -20*t*y^2;
f_y = @(t,y) -40*t*y;
t0 = 0;
T = 1;
y0 = 1;
h = 0.01;
tol = 1e-8;
N = 100;
y = imp_euler(f,f_y,t0,T,y0,h,tol,N)
plot(t0:h:T,y)
hold on
plot(t0:h:T,1./(1+10*(t0:h:T).^2))
function y = imp_euler(f,f_y,t0,T,y0,h,tol,N)
t = t0:h:T;
n = length(t);
y = zeros(n,1);
y(1) = y0;
for k = 1:n-1
g = @(z) z - y(k) - h*f(t(k+1),z);
gp = @(z) 1 - h*f_y(t(k+1),z) ;
disp('Vor Newton call.')
y(k+1) = newton(g,gp,y(k),tol,N);
disp('nach Newton Call.')
end
end
function sol=newton(f,fp,x0,tol,N)
i=0;
sol = zeros(N,2);
fc=abs(f(x0));
while fc > tol
xc = x0 - (f(x0)/fp(x0));
fc=abs(f(xc));
x0 = xc;
i=i+1;
if (i>N)
fprintf('Method failed after %d iterations. \n',N);
break
end
end
sol = x0;
end

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by