help fixing the Heuns method

4 views (last 30 days)
Sonja Cortes
Sonja Cortes on 27 Mar 2020
Edited: Jim Riggs on 27 Mar 2020
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

Answers (1)

Jim Riggs
Jim Riggs on 27 Mar 2020
Edited: Jim Riggs on 27 Mar 2020
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 Comments
Jim Riggs
Jim Riggs on 27 Mar 2020
Edited: Jim Riggs on 27 Mar 2020
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 on 27 Mar 2020
Edited: Jim Riggs on 27 Mar 2020
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.

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!