Euler VS Runge kutta – Circle – Matlab code

In the context of this exercise, lets suppose having the two dependent following function:

circle-matlab-code

With initial conditions

circle-matlab-code

The equation of a circle look like the following

circle-matlab-code

where the r is a constant.

If we add our equations at 0, the following will be true

circle-matlab-code

We are now going to be using Euler methods and the Runge-Kutta to model a circle of a radius 1 using the above mentioned equations.

Forward Euler method

Based on the equations

circle-matlab-code

We can write the following based on our example.

circle-matlab-code

Here is the Matlab code deriving from the above.

clc;

clear all;

h=0.2; % value of the step, you can vary this to see how it affect the graph

N=51;

x=zeros(1,N);

y=zeros(1,N);

x(1)=1; %initial value x(0)

y(1)=0; %initial value y(0)

%iteration with 50 step

% Forward Euler method

for i=1:N

x(i+1)=x(i)+h*(y(i));

y(i+1)=y(i)-h*(x(i));

end

angle=linspace(0,2*pi); % used to plot the unit circle

plot(x,y,'r',sin(angle),cos(angle),'b')

And the following plot results from it.

10

Backward Euler method

Just like the previous method we will start from the formula, then the derivation according to this exercises, the Matlab code and then the graph

Formula

circle-matlab-code

Derivation

circle-matlab-code

we need to arrange the above equation for it to look like the following (substituting and putting n+1 elements on one side and n elements on the other side )

circle-matlab-code

Matlab code

clc;

clear all;

h=0.2;

N=51;

x=zeros(1,N);

y=zeros(1,N);

x(1)=1; %initial value x(0)

y(1)=0; %initial value y(0)

%iteration with 50 step

% Forward Euler method

for i=1:N

x(i+1)=(x(i)+h*(y(i)))/(1+h^2);

y(i+1)=(y(i)-h*(x(i)))/(1+h^2);

end

angle=linspace(0,2*pi); % used to plot the unit circle

plot(x,y,'r',sin(angle),cos(angle),'b')

legend('Euler', 'Unit circle')

Plot

10

Semi-implicit Euler method

Derivation

circle-matlab-code

Matlab code

clc;

clear all;

h=0.2;

N=51;

x=zeros(1,N);

y=zeros(1,N);

x(1)=1; %initial value x(0)

y(1)=0; %initial value y(0)

%iteration with 50 step

% Forward Euler method

for i=1:N

x(i+1)=x(i)+h*(y(i));

y(i+1)=y(i)-h*(x(i+1));

end

angle=linspace(0,2*pi); % used to plot the unit circle

plot(x,y,'r',sin(angle),cos(angle),'b')

legend('Euler', 'Unit circle')

Plot

10

Runge Kutta method

Derivation

runge-kutta

Working the same way with y

runge-kutta

Matlab code

 clc;
clear all;
h=0.2;
N=51;
x=zeros(1,N);
y=zeros(1,N);
x(1)=1; %initial value x(0)
y(1)=0; %initial value y(0)
%iteration with 50 step
% Forward Euler method
for i=1:N
xk1=y(i);
yk1=-x(i);
xk2=y(i)+0.5*h*yk1;
yk2=-(x(i)+0.5*h*xk1);
x(i+1)=x(i)+h*xk2;
y(i+1)=y(i)+h*yk2;
end
angle=linspace(0,2*pi); % used to plot the unit circle
plot(x,y,'r',sin(angle),cos(angle),'b')
legend('Range Kutta', 'Unit circle')

Plot

runge-kutta-matlabcode

cad exercises
  • Cameron Dowd

    Can you please explain how you arrived at the various equations for the Runge-Kutta method?