Heun's Method not graphing right

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')

[Walter Roberson]

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.

Answer 1

[Anish Mitra]

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