Why does my graph look funny?
1 view (last 30 days)
Show older comments
Whenever I try to graph this, something tells me that there's something wrong. Shouldn't Euler's method, Heun's method, etc graph look similar to one another? Here's what I have:
f=inline('t.^2*cos(y)','t','y'); %Define the function
t=linspace(0,10,100); %Mesh points
y0=0.1; %Initial condition
ye=EulerODE(f,t, y0) %Solve using Euler's method
yem=EulerModODE(f, t, y0); %Solve using Euler's method
yh=HeunODE(f, t, y0); %Solve using Heun's method
yrk4=RK4(f, t, y0); %Solve using Runge-Kutta Method
yam=AdamsM(f, t, y0); %Solve using Adams-Moulton Method
plot(t, ye, 'Linewidth',2)
hold on
plot(t, yem, 'co', 'Linewidth', 2)
hold on
plot(t, yh, 'bx', 'Linewidth', 2)
hold on
plot(t, yrk4, 'k+', 'Linewidth', 2)
hold on
plot(t, yam, 'm*', 'Linewidth', 2)
hold off
Function EulerODE, EulerModODE, HeunODE etc. were taken from a MatLab CD provided from my text book. Below are some of the functions to give you an idea:
function x = EulerODE(f, t, x0)
%
% EULERODE Uses Euler's method to solve a first-order ODE given in the form
% xdot = f(t, x) subject to initial condition x0.
%
% x = EulerODE(f, t, x0) where
%
% f is an inline function representing xdot = f(t, x).
% t is a vector representing the mesh points.
% x0 is a scalar representing the initial value of x.
%
% x is the vector of x values at the mesh points in t.
x = 0*t; % Preallocate x
x(1) = x0;
for n = 1:length(t)-1
h = t(n+1) - t(n);
x(n + 1) = x(n) + h*f(t(n), x(n));
end
function x = EulerModODE(f, t, x0)
%
% EULERMODODE Uses the modified Euler's method to solve a first-order ODE
% given in the form xdot = f(t, x) subject to initial condition x0.
%
% x = EulerModODE(f, t, x0) where
%
% f is an inline function representing xdot = f(t, x).
% t is a vector representing the mesh points.
% x0 is a scalar representing the initial value of x.
%
% x is the vector of x values at the mesh points in t.
x = 0*t; %Preallocate x
x(1) = x0;
for n = 1:length(t)-1
h = t(n+1) - t(n);
x(n + 1) = x(n) + h*f(t(n) + h/2, x(n)+ h/2*f(t(n), x(n)));
end
Any help would be appreciated. Thanks!
0 Comments
Answers (1)
Grzegorz Knor
on 3 Dec 2011
Compare your results with MATLAB ODE functions, e.g.:
ode45(f,[0 10],0.1)
If you increase the number of mesh points:
t=linspace(0,10,1000); %Mesh points
results will be comparable.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!