Euler's method for two first order differential equations?

I was trying to solve two first order differential equations like below using the Euler's method and plot two graphs with x and y as a function of t.

The differential equations are:

dxdt = @(x,t) -1.*y-0.1.*x; 
dydt = @(y,t) x-0.5.*y;

I tried this script below:

a=0; %initial time
b=2000; %final time

h = 0.01; % time step
N = (b-a)./h;
t=a:h:b;
t(1)=a;
x(1)=0;
y(1)=0;

dxdt = @(x,t) -1.*y-0.1.*x; 
dydt = @(y,t) x-0.5.*y;
for n=1:N
    
    x(n+1)=x(n)+h.*dxdt(x(n),t(n));
    y(n+1)=y(n)+h.*dydt(y(n),t(n)); 
    t(n+1)=a+n*h;
end
plot(t,x,'-',t,y,'--')

But there was an error saying

Unable to perform assignment because the left and right sides have a different
number of elements.

Error in Q6 (line 15)
    x(n+1)=x(n)+h.*dxdt(x(n),t(n));

I did not quite understand why.

Could anyone help me with this please?

Thanks in advance.

[darova]

Where is t?

[Yuhong Jin]

t is in the differential equations, dxdt and dydt.

I only have the range of t and relationships between dxdt, dydt and x, y.

Answer 1

[Jim Riggs]

Try preallocating x and y

nsteps = (b-a)/h + 1;  % this is the number of elements in t(a:h:b) (It is = N+1)
x=zeros(1,nsteps);
y=zeros(1,nsteps);

Also, when you define t as

t=a:h:b;

This creates a vector with T(1)=a and t(end)=b, so the following line

t(1)=a

is not necessary.

Inside your for loop, t has already been defined, above, as t=a:h:b, so you don't need to re-define it.

x and y need to be preallocated to size "nsteps", which is 1 greater than N, in order for your loop to work.

[Jim Riggs]

I just noticed that your anonymous functions do not have the propper form.

You need to supply both functions with both x and y as inputs;

dxdt = @(x,y) -1*y-0.1*x;
dydt = @(x,y) x-0.5*y;

Then inside your loop:

x(n+1) = x(n) + h*dxdt(x(n), y(n));
y(n+1) = y(n) + h*dydt(x(n), y(n));

[Chris Horne]

I'm trying to run a similar code and get the same error: "Unable to perform assignment because the left and right sides have a different number of elements. The line of code in bold

Error in EulerHorne (line 99)

x(n+1) = x(n) + h*dxdt(x(n), y(n));

a=0

b=2000

h=0.01; % step size

nsteps = (b-a)./h + 1; % this is the number of elements in t(a:h:b) (It is = N+1)

t=a:h:b;

x(1)=0.1; y(1)=0.1; % Initial x,y

x=zeros(1,nsteps); % preallocate x and y

y=zeros(1,nsteps); % preallocate x and y

dxdt = @(x,y) x./(2*(t+1))-2*t*y; % RHS of x equation

dydt = @(x,y) y./(2*(t+1))-2*t*x; % RHS of x equation

for n=1:nsteps

x(n+1) = x(n) + h*dxdt(x(n), y(n));

y(n+1) = y(n) + h*dydt(x(n), y(n));

end

figure;

plot(x,y) % Plot x vs y

xlabel('x')

ylabel('y')

load gong

sound(y,Fs)

[Torsten]

Why did you return to the notation of integrating all equations separately ?

Say you have 100 equations to solve. Do you want to write the line

dxdt = @(x,y) x./(2*(t+1))-2*t*y; % RHS of x equation

and the line

x(n+1) = x(n) + h*dxdt(x(n), y(n));

100 times ?

The code you wrote here

https://de.mathworks.com/matlabcentral/answers/1678849-how-to-specify-a-3-element-column-vector-in-euler-s-method-for-ode?s_tid=srchtitle

will work perfectly in the above case.

[Chris Horne]

%Tristen, OK. You're right. Back to the 3D code that worked. 
%I'm getting ARRAY SIZE error for the 2D version INSIDE THE LOOP
% Can you help?

%Tristen, OK. Youre right. Back to the 3D code that worked. I'm getting error for the 2D version INSIDE THE LOOP
% Can you help?
% solve u'(t) = F(t,u(t)) where u(t)= u1./(2*(t+1))-2*t*u2, u2./(2*(t+1))-2*t*u1;
% Euler's Method
% Initial conditions and setup
neqn = 2; % set a number of equations variable
h=input('Enter the step size for 2D system:  ') % step size will effect solution size
t=(0:h:4).';%(starting time value 0):h step size
nt = size(t,1); % size of time array
%(the ending value of t2 );  % the range of t
% define the function vector, F
F = @(t,u)[u1./(2*(t+1))-2*t*u2,u2./(2*(t+1))-2*t*u1]; % define the function 'handle', F
% with hard coded vector functions of time
u = zeros(nt,neqn); % initialize the u vector with zeros
w=input('Enter the intial vector values of 2 components using brackets [u(0)= 1,v(0)=0]:  ')
u(1,:)= w; % the initial u value amd the first column
%n=numel(u); % the number of u values
% The loop to solve the ODE (Euler's Method)
for i = 1:nt-1
  u(i+1,:) = u(i,:) + h*F(t(i),u(i,:)); % Euler's formula for a 2D vector function F
end
figure;
plot(u1,u2) % Plot u vs v 
xlabel('u')
ylabel('v')

[Chris Horne]

Sorry Torsten not Tristen

[Torsten]

neqn = 2; % set a number of equations variable
h=0.01;%input('Enter the step size for 2D system:  ') % step size will effect solution size
t=(0:h:4).';%(starting time value 0):h step size
nt = size(t,1); % size of time array
%(the ending value of t2 );  % the range of t
% define the function vector, F
F = @(t,u)[u(1)./(2*(t+1))-2*t*u(2),u(2)./(2*(t+1))-2*t*u(1)]; % define the function 'handle', F
% with hard coded vector functions of time
u = zeros(nt,neqn); % initialize the u vector with zeros
w=[1 2];%input('Enter the intial vector values of 2 components using brackets [u(0)= 1,v(0)=0]:  ')
u(1,:)= w; % the initial u value amd the first column
%n=numel(u); % the number of u values
% The loop to solve the ODE (Euler's Method)
for i = 1:nt-1
    u(i+1,:) = u(i,:) + h*F(t(i),u(i,:)); % Euler's formula for a 2D vector function F
end
figure;
plot(u(:,1),u(:,2)) % Plot u vs v 
%plot(t,u)
xlabel('u')
ylabel('v')

[Chris Horne]

Sorry I had an error in my two functions.

% solve u'(t) = F(t,u(t)) where u(t)= sqrt(t+1)*cos^2(t),sqrt(t+1)*sin^2(t) ;
% Euler's Method
% Initial conditions and setup
neqn = 2; % set a number of equations variable
h=input('Enter the step size for 2D system:  ') % step size will effect solution size
t=(0:h:4).';%(starting time value 0):h step size
nt = size(t,1); % size of time array
%(the ending value of t2 );  % the range of t
% define the function vector, F
F = @(t,u)[sqrt(t+1)*cos.^2(t),sqrt(t+1)*sin.^2(t)]; % define the function 'handle', F
% with hard coded vector functions of time
u = zeros(nt,neqn); % initialize the u vector with zeros
w=input('Enter the intial vector values of 2 components using brackets [u(0)= 1,v(0)=0]:  ')
u(1,:)= w; % the initial u value amd the first column
%n=numel(u); % the number of u values
% The loop to solve the ODE (Euler's Method)
for i = 1:nt-1
  u(i+1,:) = u(i,:) + h*F(t(i),u(i,:)); % Euler's formula for a vector function F
end
%fprintf('="U"\n\t   %0.01f',u);

[Torsten]

cos.^2(t) and sin.^2(t) are wrong.

cos(t).^2 and sin(t).^2 are correct (although you don't need the elementwise multiplication (.) here since the values given to F are always scalars).

[Chris Horne]

I figured it out. The function definition delimiters typo. The code produces then 2D plot of the approximating functions.

[Chris Horne]

Is the term 'forward Euler' the same as 'Euler' in terms of the algorithm? Here is my method for solving 3 equaitons as a vector:

% This code solves u'(t) = F(t,u(t)) where u(t)= t, cos(t), sin(t)

% using the FORWARD EULER METHOD

% Initial conditions and setup

neqn = 3; % set a number of equations variable

h=input('Enter the step size: ') % step size will effect solution size

t=(0:h:4).';%(starting time value 0):h step size

nt = size(t,1); % size of time array

%(the ending value of t ); % the range of t

% define the function vector, F

F = @(t,u)[t,cos(t),sin(t)]; % define the function 'handle', F

% with hard coded vector functions of time

u = zeros(nt,neqn); % initialize the u vector with zeros

v=input('Enter the intial vector values of 3 components using brackets [u1(0),u2(0),u3(0)]: ')

u(1,:)= v; % the initial u value and the first column

%n=numel(u); % the number of u values

% The loop to solve the ODE (Forward Euler Algorithm)

for i = 1:nt-1

u(i+1,:) = u(i,:) + h*F(t(i),u(i,:)); % Euler's formula for a vector function F

end

Answer 2

[Eng]

% Define the differential equation dydx = @(x,y) x + y;

% Define the initial condition y0 = 1;

% Define the step size and interval h = 0.1; x = 0:h:1;

% Initialize the solution vector y = zeros(size(x)); y(1) = y0;

% Apply the Modified Euler Method for i = 1:length(x)-1 k1 = dydx(x(i), y(i)); k2 = dydx(x(i+1), y(i)+hk1); y(i+1) = y(i) + h/2(k1+k2); end

% Plot the solution plot(x,y) xlabel('x') ylabel('y')