Any ideas why Euler method isn't running?

1 Ansicht (letzte 30 Tage)
Kaylene Widdoes
Kaylene Widdoes am 17 Feb. 2016
Beantwortet: Geoff Hayes am 18 Feb. 2016
I have the following equation: (x^3)y' + 20*(x^2)y = x with y(2)=0
I have to get the Euler code to run and execute. So far this is my code, but the graph looks wrong, and the approximations aren't even showing up.
Script file:
function[T,Y] = euler(f,a,b,ya,h)
a = 0; b = 1;
T = a:h:b;
Y = zeros(1,length(T));
Y(2) = ya;
for k = 1 : length(T)-1
Y(k+1) = Y(k) + h*f(T(k),Y(k));
end
Code:
%%Part 1 - Euler
% y'(t) = 3*y(t)
% A function file euler.m is needed; you'll do this in the homework
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]=euler(f,a,b,ya,h);
% fine solution
h = .001;
ya = 0;
[T2,Y2]=euler(f,a,b,ya,h);
% finer solution
h = .0001;
ya = 0;
[T3,Y3]=euler(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.001')
title('The Euler Method with 3 meshes')

Antworten (1)

Geoff Hayes
Geoff Hayes am 18 Feb. 2016
Kaylene - if I run your code, I notice that your Y array output from euler are all NaNs. Looking closely at this function
function[T,Y] = euler(f,a,b,ya,h)
a = 0; b = 1;
T = a:h:b;
Y = zeros(1,length(T));
Y(2) = ya;
for k = 1 : length(T)-1
Y(k+1) = Y(k) + h*f(T(k),Y(k));
end
the a and b inputs are overwritten with 0 and 1 respectively. Why? If I remove these initializations and change the initial value for Y (so Y(1) and not Y(2)) then the code becomes
function[T,Y] = euler(f,a,b,ya,h)
T = a:h:b;
Y = zeros(1,length(T));
Y(1) = ya;
for k = 1 : length(T)-1
Y(k+1) = Y(k) + h*f(T(k),Y(k));
end
which (at least) initializes Y to be something that can be drawn to your figure. The output using the different step sizes appears to be near identical which may not be what you are trying to show.

Kategorien

Mehr zu Linear Algebra finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by