Function as an integrator

Hi everybody,

I'm computing a numerical solution by Euler's method for a second order differential equation, I have defined the function for the differential equation but I also need to define an additional one as an integrator (like ode45), but by Euler's method. I've tried some computations and agorithms but it does not work. Please, could you help me?

This is my code:

% Sampling interval
h = 0.1;  
%h = 0.05;
%h = 0.01;

t0 = 0; % Initial time
tf = 4; % Final time

t  = [t0:h:tf];  % Time vector

% Initial value problem
th0 = pi/40;  % theta(t0)
om0 = 0;      % ommega(t0)
yu = [th0 om0];  

e = euler(@yprime,t,h,yu); % Want to call euler's function as an integrator but does not work and has an issue with * multiplication

function [yp] = yprime(t,y)
L = 0.61;
G = 9.81;
yp = zeros(2,1);

yp(1) = y(2); % Angle theta [rad]
yp(2) = (-G/L)*sin(y(1)); % Angular velocity omega [rad/s]

end

% and my Euler's function as a script:

function eu = euler(yp,t,h,y)
% Function that computes the values of the angle theta followinf the Euler
% method

eul = y + yp*h;

eu = eul;

Answer 1

[Walter Roberson]

eul = y + yp(t,y)*h;

[Hugo Hernández Hernández]

Thank you so much, that solved the problem with multiplication in euler's function, but as an output I'm getting just a 2x2 matrix for e, and should be 2x41 the lenght of t. I tried asgning a zeros matrix in euler's function with t length but didn't work, should this zeros matrix be implemented in euler's or differential equation function? The differential eq. function worked a expected with an ode45 integrator, but I do not know how to get this set of values using this euler's function.

[Walter Roberson]

eul = y(:) + yp(t,y)*h;

Then you will get a proper 2 x 1 output.

2 x length(t) is not to be expected, as your code is independent of t. You do not loop over t values, and you ignore t in your calculations.

[Hugo Hernández Hernández]

Ok, many thanks.

[Hugo Hernández Hernández]

Hi again Walter, I'm still having problems with that, for the euler's function I have defined as follows:

function eu = euler(yp,t,h,y)
% Function that computes the values of the angle theta following the Euler
% method

for i = 1:length(t)-1
    eu1(i+1) = y(i)+h*yp(t(i));
%eul = y(:) + yp(t,y)*h;
end

eu = eul;

And then trying to implement in the main code but I get errors like:

Not enough input arguments.

Error in test>yprime (line 34)

yp(1) = y(2); % Angle theta [rad]

Error in euler (line 6)

eu1(i+1) = y(i)+h*yp(t(i));

Error in test (line 27)

e = euler(@yprime,t,h,yu);

I've seen many examples with euler's method but seem that any is calling a function like the defined in my main script, what could be the problem?

Thanks in advance.

[Walter Roberson]

    eu1(i+1) = y(i)+h*yp(t(i));

That line of code does not appear in anything you have posted so far. It is failing to pass y values into yp. It is also trying to work with one y value at a time, which is incorrect as yp needs both y values and returns two outputs.

[Hugo Hernández Hernández]

Got it, thanks!