Solving ODE using Euler Method and 4th Order Runge Kutta Method

I am trying to learn how to solve differential equations provided the intial conditions, I have already made the matlab code for both the euler and runge kutta method as follows:

%Euler method

function y=elrl(t,x,n,h)
%t:Time
%x:Variables
%n:Number of variables
%h:Step
d=f(t,x);
for i = 1:n
    p(i)=x(i)+h*d(i);   
end
d=f(t+h,p);
for i = 1:n
    q(i)=x(i)+h*d(i); 
    y(i)=0.5*(p(i)+q(i));
end

%Runge Kutta method

function y=rktl(t,x,n,h)
%t:Time
%x:Variables
%n:Number of variables
%h:Step
   
k0=f(t,x);
k1=f(t+0.5*h,x+k0*0.5*h);
k2=f(t+0.5*h,x+k1*0.5*h);
k3=f(t+h,x+k2*h);

y=x+h/6*(k0+2*k1+2*k2+k3);

My problem is that I don't know how to define the function using the initial conditions, I have the following system, can anybody help ?

[Torsten]

What are the x-variables in your code ?

[Adnane Ait Zidane]

Hello, thank you for replying.

The x variables are just the y functions used to describe the euler method and runge kutta method.

[James Tursa]

None of the methods is coded correctly. In the first place, you have three 1st order differential equations, so you will need a 3-element state vector at each time point to hold the solution. But you have coded a 1-element scalar as your solution at each time step. Secondly, your Euler schemes apparently attempt to use the same incoming state for all steps when they should be using a calculated state iteravely at each step. And it looks like your RK4 scheme only takes one step. Can you show us the code you have written for your f( ) function? I suspect that has errors also.

[Adnane Ait Zidane]

I thought so aswell, I have done the following to use the Euler method. The problem is that I do not have an explicit expression for the function I want to appriximate. The only information I have is in the picture linked above.

New code:

y0=-1;y1=0;y2=1;

t=0;
n=1000;
dt=0.01;

for i=1:n
    y0=y0+dt*y1;
    y1=y1+dt*(-y0);
    y2=y2+dt*(-y2);
    
    t=t+dt;
    p(i,:)=[y0 y1 y2];
end
figure
plot(p)

[Torsten]

Yes, that's the explicit Euler method for your system.

But you should replace the plot command by

plot (transpose(0.01:0.01:10),p)

Do you see where the initial conditions are set ?

[Adnane Ait Zidane]

Yes, from the figure it's clear where y0 y1 and y2 are starting from. In the y axis we can read the values, am i right?

I see why I should use the other command plot, it makes for reading values easier, thank you.

[James Tursa]

These two lines are not quite correct:

y0=y0+dt*y1;
y1=y1+dt*(-y0);

The problem is that you are using the current y1 to propagate the y0 variable (which is correct), but you are using the next y0 to propagate the y1 variable (which is incorrect). You should do something like this instead:

y0save = y0;
y0=y0+dt*y1;
y1=y1+dt*(-y0save);

That way you are using the current states on the right hand side consistently.

[Adnane Ait Zidane]

How is possible to define the function in a different script and use the euler method in other another and combine them? I have just started using matlab so i am not very proefficient.

Thank you.

[Adnane Ait Zidane]

Am I understanding correctly?

I have to save the initial value of y0 to use it to calculate y1 ?

Answer 1

[James Tursa]

You pretty much have the Euler scheme worked out, so I will help you with a vector formulation. Take this code:

for i=1:n
    y0=y0+dt*y1;
    y1=y1+dt*(-y0); % but as noted this line isn't quite correct
    y2=y2+dt*(-y2);
    t=t+dt;
    p(i,:)=[y0 y1 y2];
end

Here you are filling p(i.:) at each step with the calculated y0, y1, and y2 variables which you calculate individually. However, you could code this directly as a vector like this (note that I have switched the indexes so that the state is a column vector):

for i=1:n-1
    p(:,i+1) = p(:,i) + dt * f(t,p(:,i));
    t=t+dt;
end

In this case the p(:,i) is a 3-element state vector of [y0;y1;y2] at the current time t, and p(:,i+1) is the 3-element state vector at the next time t+dt. It follows your scheme where p(1,i) is y0, p(2,i) is y1, and p(3,i) is y2. The only thing missing is the vector derivative function f, which could be simply this:

f = @(t,y) [y(2);-y(1);-y(3)]; % where y(1)=y0, y(2)=y1, y(3)=y2

Why did I switch the indexing? By having the states in columns, your derivative function will match what the MATLAB supplied ode functions such as ode45 expect, and it will be easy for you to double check your results by calling ode45 using the same f function. Also, it will be easier to take this vector formulation and extend it to the Modified Euler method and the RK4 scheme.

[Adnane Ait Zidane]

I have done the following, I get an error though. Both are 1x1 doubles so where is the problem coming from?

%Euler Method

function y=elrl(y,h)
n=1000;                     %Number of iterations
t=0;                        %Initial value
dt=0.01;                    %Step size
y0=-1;y1=0;y2=1;

for i=1:n
    y0save=y0;
    y0=y0+dt*y1;
    y1=y1+dt*(-y0save);
    y2=y2+dt*(-y2)
    p(i,:)=[y0 y1 y2];
end

for i=1:n-1
    f = @(t,y) [y(2);-y(1);-y(3)];
    p(:,i+1)=p(:,i)+dt*f(t,p(:,i));
    t=t+dt;
end

figure
plot (transpose(0.01:0.01:10),p)

Answer 2

[Torsten]

function main
  tstart = 0.0;
  tend = 10.0;
  h = 0.01;
  T = (tstart:h:tend).';
  Y0 = [-1 0 1]; 
  f = @(t,y) [y(2) -y(1) -y(3)];
  Y = euler_explicit(f,T,Y0);        % Euler explicit solution
  plot(T,Y)
  hold on
  Y = runge_kutta_RK4(f,T,Y0);
  plot(T,Y)                           % Runge Kutta solution
  %hold on
  %plot(T,[-cos(T) sin(T) exp(-T)])  % Analytical solution
end

function Y = euler_explicit(f,T,Y0)
  N = numel(T);
  n = numel(Y0);
  Y = zeros(N,n);
  Y(1,:) = Y0;
  for i = 2:N
    t = T(i-1);
    y = Y(i-1,:);
    h = T(i) - T(i-1);
    Y(i,:) = y + h*f(t,y);
  end
end

function Y = runge_kutta_RK4(f,T,Y0)
  N = numel(T);
  n = numel(Y0);
  Y = zeros(N,n);
  Y(1,:) = Y0;
  for i = 2:N
    t = T(i-1);
    y = Y(i-1,:);
    h  = T(i) - T(i-1);
    k0 = f(t,y);
    k1 = f(t+0.5*h,y+k0*0.5*h);
    k2 = f(t+0.5*h,y+k1*0.5*h);
    k3 = f(t+h,y+k2*h);
    Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
  end
end 

You must improve your programming skills.

So I hope you do not only copy the code and submit it as your homework, but also try to understand what's going on here.

[Adnane Ait Zidane]

I would like to understand the logic behind putting the expression of the function into that order, it would be helpful if you can explain.

Thank you.

[Torsten]

The main principle you have to follow is that the variables you use in a certain line of a function must be defined or calculated in the lines before.

So you should first consider what you want the function to supply as output, what variables are known (input to the function) and what variables you need to calculate from these known quantities to produce the desired output.