MATLAB Answers

help fixing the Heuns method

5 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

  0 Comments

Sign in to comment.

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
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
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.

Tags


Translated by