Trying to make euler's method program work

1 view (last 30 days)
ce
ce on 5 Dec 2022
Edited: Torsten on 6 Dec 2022
function [x,y] = euler_for(x0,y0,x_end,h,fcn)
%Solve the ivp
% y' = f(x,y), x0 <= x <= b, y(x0) = y0
%Use euler's method with a stepsize h. The
%user must supply a program with some name, say, deriv,
%and a first line of the form
% function ans = deriv(x,y)
%A sample call would be
% [t,z] = euler_for(t0,z0,b,delta,'deriv')
%
%output:
%the routine eulercls will return two vectors, x and y.
%The vector x will contain the node points
% x(1) = x0, x(j) = x0 + (j-1)*h, j=1,2,...,N
%with
% x(N) <= x_end - h, x(N) + h > x_end - h
%The vector y will contain the estimates of the solution Y
%at the node points in x.
%
n = fix((x_end-x0)/h)+1;
x = linspace(x0,x_end,n);
y = zeros(n,1);
y(1) = y0;
for i = 2:n
y(i) = y(i-1)+h*feval(fcn,x(i-1),y(i-1));
end
I'm not sure what is is meant by the function that I need to supply. I know what function I'm supposed to use, but I don't know how to add it to the program. I tried making a function called fcn but when I run the program it says there is an error in fcn, although it doesn't say what it is. This is my fcn code:
function f = fcn(x,y)
f = (cos(atan(x)))^2;
end

Answers (1)

Alan Stevens
Alan Stevens on 5 Dec 2022
Like this:
% Arbitrary values here - insert your own
x0 = 0;
x_end = 10;
h = 0.1;
y0 = 1;
n = fix((x_end-x0)/h)+1;
x = linspace(x0,x_end,n);
y = zeros(n,1);
y(1) = y0;
for i = 2:n
y(i) = y(i-1)+h*fcn(x(i-1));
end
plot(x,y),grid
xlabel('x'),ylabel('y')
function f = fcn(x) % Your function doesn't depend on y
f = (cos(atan(x)))^2;
end
  2 Comments
Torsten
Torsten on 6 Dec 2022
Edited: Torsten on 6 Dec 2022
x0 = 0;
x_end = 10; % Integration interval [x0,x_end]
h = 0.1; % Step size of integration
fcn = @(x,y) (cos(y))^2; % Y' = cos(Y)^2
y0 = 0; % Y(x0) = 0
[x,y] = euler_for(x0,y0,x_end,h,fcn);
hold on
plot(x,y),grid
xlabel('x'),ylabel('y')
syms y(x)
Ysol(x) = dsolve(diff(y,x)==cos(y)^2,y(0)==0);
fplot(Ysol,[x0 x_end])
hold off
function [x,y] = euler_for(x0,y0,x_end,h,fcn)
n = fix((x_end-x0)/h)+1;
x = linspace(x0,x_end,n);
y = zeros(n,1);
y(1) = y0;
for i = 2:n
y(i) = y(i-1)+h*fcn(x(i-1),y(i-1));
end
end

Sign in to comment.

Categories

Find more on Mathematics in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!