Heun's Method not graphing right
조회 수: 2 (최근 30일)
이전 댓글 표시
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
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
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
댓글 수: 0
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!