MATLAB Answers

Heun's Method not graphing right

3 views (last 30 days)
Kaylene Widdoes
Kaylene Widdoes on 18 Feb 2016
Edited: Anish Mitra on 22 Feb 2016
I have the following equation: (x^3)y' + 20*(x^2)y = x with y(2)=0 that I'm trying to graph from 2,10
The graph is running really weird. Any ideas?
Script:
function [T,Y] = heun(f,a,b,ya,h)
m=(b-a)/h;
T = zeros(1,m+1);
Y = zeros(1,m+1);
T(1) = a;
Y(1) = ya;
for j=1:m,
k1 = feval(f,T(j),Y(j));
p = Y(j) + h*k1;
T(j+1) = a + h*j;
k2 = feval(f,T(j+1),p);
Y(j+1) = Y(j) + h*(k1 + k2)/2;
end
end
Code:
%%Part 2 - Heun's Method
clear all
% define the problem: function f and domain
f = @(t,y) (1/(t^2)) - ((20*y)/t);
a = 2; b = 10;
% exact solution, using a fine grid
t = a:.0001:b;
y = (1./(19.*t)) - (524288./(19.*(t.^20))); % this is a vector of values, not a function
% coarse solution
h = .01;
ya = 0;
[T1,Y1]=heun(f,a,b,ya,h);
% fine solution
h = .001;
ya = 0;
[T2,Y2]=heun(f,a,b,ya,h);
% finer solution
h = .0001;
ya = 0;
[T3,Y3]=heun(f,a,b,ya,h);
plot(t,y,'k',T1,Y1,'bo-',T2,Y2,'ro-',T3,Y3,'go-')
legend('Exact','h=0.01','h=0.001','h=0.0001')
title('The Heun Method with 3 meshes')

  1 Comment

Walter Roberson
Walter Roberson on 18 Feb 2016
Please document your code, or at least use descriptive variable names. I cannot tell if the difficulties that I see are intentional or not because nothing is clear about what your code is intended to do or about how it thinks it is doing it.

Sign in to comment.

Answers (1)

Anish Mitra
Anish Mitra on 22 Feb 2016
Edited: Anish Mitra on 22 Feb 2016
Hi,
I am assuming that when you say "The graph is running really weird", it is referring to the fact that the graph has a thick green line and not to the actual data values shown by the plot.
I believe that 'y' represents the closed-form solution of the equation, and that you want to compare the numerical solution (with different step times - in terms of step size of 'h').
Since Y3 is very finely spaced, the markers dominates the graph and it is difficult to observe the other lines. If you plot without the markers, you can see that the lines are close together.
>> plot(t,y,'k',T1,Y1,'b-',T2,Y2,'r-',T3,Y3,'g-')
Also, to check the error in each case, you can plot the difference between the closed-form solution and the numerical solution.
1. To compare y and Y3
>> figure; plot(t,y - Y3);
2. To compare y and Y1 (or Y2), we will need to first resample y to the size of Y1/Y2.
>> y1 = interp1(t,y,T1);
>> figure; plot(T1,y1 - Y1);
Hope this helps!
Regards,
Anish

  0 Comments

Sign in to comment.

Tags