2D - EDO: Numerical Methods

The different numerical methods for approximating the solution of an ODE always assume that it is a first order one When our problem is defined by an ODE of bigger order we have to transform it to a system of ODEs of order one. On this process you have to introduce as many new variables as whe n is the order of the ODE. For exemple, if we have a second order ODE:

Then, the ODE is writen in vectorial form like where is a 2-dim variable. Thus, the original variables .

% Define the ODE

f=@(t,X) [X(2); -X(1)]; %derivatives funtion

% initial values

X0=[1;0];

tini=0;

tfin=15;

% Call the numerical methods to compare

N=500;

%dt=(tfin-tini)/N;

[te,Xe]=euler(f,[tini,tfin],X0,N);

[tr,Xr]=RK4(f,[tini,tfin],X0,N);

% opcions=odeset('RelTol',1.e-12,'AbsTol',1.e-12);

% [t,X]=ode45(f,[tini tfin],X0,opcions);

t=tr; %choose the method to plot

X=Xr;

% Plot the solution

% Characteristics of the plot window

fig1=figure;

set(fig1,'Name','INTEGRACIO EDOs',...

'units', 'pixels','Position', [200, 50, 400, 600],...

'NumberTitle','off','Resize','on');

% --------------(t,x(t)) the x's as a function of time

subplot(2,2,1);

plot(t,X(:,1));

xlabel('temps');

ylabel('x');

% --------------(t,y(t)) the y's as a function of time

subplot(2,2,2);

plot(t,X(:,2));

xlabel('temps');

ylabel('y');

% --------------(x(t),y(t)) Phase portrait

subplot(2,2,3);

plot(Xe(:,1),Xe(:,2));

xlabel('x-Euler');

ylabel('y-Euler');

axis equal;

subplot(2,2,4);

plot(Xr(:,1),Xr(:,2));

xlabel('x-RK4');

ylabel('y-RK4');

axis equal;