help fixing the Heuns method

4 ビュー (過去 30 日間)
Sonja Cortes
Sonja Cortes 2020 年 3 月 27 日
編集済み: Jim Riggs 2020 年 3 月 27 日
So I need some help with this heuns method, There is something wrong in my code but I don't know where, since when i plot it the Euler method i have is more close to the function than the Heun
f: function, i.e. right side of equation y '= f (x, y)
t0 and y0: initial values such that y (x0) = y0
no_steg: number of Heun steps The Heun function should calculate
lengde: number that indicates the length of an HEun step
t0 = 0; %start value of x
y0 = 2; %start value of y
no_steg = 3;
lengde = 0.2;
t = 0:lengde:no_steg;
f= @(t,y) -2*t*y; % function
[xverd, tilnaer] = Heun_met (f, t0, y0, no_steg, lengde)
plot(xverd, tilnaer)
hold on
y1= y0;
for s=1: length(t)-1
t0(s+1) = t0(s) + lengde;
y1(s+1)= y1(s) + lengde*f(t0(s), y1(s)); %Beregner et euler steg for startverdi problemet
end
plot(xverd, y1)
y2 = 2*exp(-t.^2);
plot(xverd, y2)
hold off
function [xverd, tilnaer] = Heun_met (f, t0, y0, no_steg, lengde)
t = 0:lengde:no_steg;
y1 = y0;
for s=1: length(t)-1
t0(s+1) = t0(s) + lengde;
y1(s+1)= y1(s) + lengde*f(t0(s), y1(s)); %Beregner et euler steg for startverdi problemet
y0(s+1) = y0(s) + (lengde/2) * (f(t0(s), y0(s)) + (f(t0(s+1) + lengde, y1(s+1)))); %Beregner heun steg
end
tilnaer = y0;
xverd = t0;
end

回答 (1 件)

Jim Riggs
Jim Riggs 2020 年 3 月 27 日
編集済み: Jim Riggs 2020 年 3 月 27 日
Heun's method is otherwise known as "explicit trapezoidal method", "improved Euler's method" or "modified Euler's method".
The way you have it,
t0(s+1) = t0(s) + lengde,
%therefore,
t0(s+1) + lengde = t0(s) + 2*lengde
therefore, "lengde" is being added twice, which is incorrect
[see comment, below]
  2 件のコメント
Jim Riggs
Jim Riggs 2020 年 3 月 27 日
編集済み: Jim Riggs 2020 年 3 月 27 日
function [t, y] = Heun_met (f, t0, y0, no_steg, lengde)
t = t0:lengde:no_steg;
y = 0.*t;
y(1) = y0;
for s=1: length(t)-1
k1 = lengde * f(t(s), y(s));
k2 = lengde * f(t(s) + lengde, y(s) + k1)
y(s+1) = y(s) + (k1 + k2)/2;
end
end
Jim Riggs
Jim Riggs 2020 年 3 月 27 日
編集済み: Jim Riggs 2020 年 3 月 27 日
The other way to code it is:
...
for s=1: length(t)-1
k1 = f(t(s), y(s));
k2 = f(t(s) + lengde, y(s) + lengde*k1)
y(s+1) = y(s) + lengde * (k1 + k2)/2;
end
...
This is equivalent to the one, above.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeProgramming についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by