Heun's Method not graphing right

5 次查看（过去 30 天）
Kaylene Widdoes
Kaylene Widdoes 2016-2-18
编辑： Anish Mitra 2016-2-22
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 个评论
Walter Roberson
Walter Roberson 2016-2-18
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Anish Mitra
Anish Mitra 2016-2-22
编辑：Anish Mitra 2016-2-22
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

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by