MATLAB Answers

help fixing the Heuns method

조회 수: 1(최근 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일
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.

댓글을 달려면 로그인하십시오.

태그

Community Treasure Hunt

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

Start Hunting!

Translated by