Question

Consider the following second-order ODE,

y"+1/4y=0

with y(0) = 1 and y'(0)=0.

Transform this unique equation into a system of two 1st-order ODEs.

Solve the obtained system for t in [0,0.6] with h = 0.2 by MATLAB using Euler Method or Improved Euler Method.

Answer 1

%Matlab code for Euler and Heun's and RK4 method for 2d ODE
clear all
close all
%functions for Euler and Huen method

f=@(t,y1,y2) y2;
g=@(t,y1,y2) (-1/4)*y1;

fprintf('Here let us consider\n')
fprintf('\tdx/dt=f(t,y1,y2)=y2\n')
fprintf('\tdy/dt=g(t,y1,y2)=-(1/4)*y1\n')
%step size
h=0.2;

%Initial guess
y1_0=1; y2_0=0;

y1_euler(1)=y1_0; y2_euler(1)=y2_0;

t_euler(1)=0; t_end=0.6;
n=(t_end-t_euler(1))/h;

%iteration for y1 and y2 using Euler method
fprintf('Iterations for y1 and y2 using Euler method\n\n')
fprintf('Initial condition for y1(0)=%f and y2(0)=%f .\n',y1_0,y2_0)
fprintf('Euler iterations\n')
for i=1:n+1
    fprintf('After %d iterations \n',i)
    t_euler(i+1)=t_euler(i)+h;
    y1_euler(i+1)=y1_euler(i)+h*f(t_euler(i),y1_euler(i),y2_euler(i));
    y2_euler(i+1)=y2_euler(i)+h*g(t_euler(i),y1_euler(i),y2_euler(i));
    fprintf('\tt_euler(%d+1)=t_euler(%d)+h=%f\n',i,i,t_euler(i+1));
    fprintf('\n\ty1_euler(%d+1)=y1_euler(%d)+h*f(t_euler(%d),y1_euler(%d),y2_euler(%d))=%f\n',i,i,i,i,i,y1_euler(i+1))
    fprintf('\ty2_euler(%d+1)=y2_euler(%d)+h*g(t_euler(%d),y1_euler(%d),y2_euler(%d))=%f\n',i,i,i,i,i,y2_euler(i+1))
end
fprintf('\n\t Using Euler method at t=%f ,y(%f)=%f\n',t_euler(end),t_euler(end),y1_euler(end))
%iteration for y1 and y2 using Huen method
fprintf('\n\nIterations for x and y using Heun method\n\n')
fprintf('Initial condition for x(0)=%f and y(0)=%f .\n',y1_0,y2_0)
y1_huen(1)=y1_0; y2_huen(1)=y2_0;
t_huen(1)=0; t_end=0.6;
n=(t_end-t_huen(1))/h;

for i=1:n+1
    fprintf('After %d iterations \n',i)
    t_huen(i+1)=t_huen(i)+h;
    k1=h*f(t_huen(i),y1_huen(i),y2_huen(i));
    l1=h*g(t_huen(i),y1_huen(i),y2_huen(i));
    fprintf('\tt_heun(%d+1)=t_heun(%d)+h=%f\n',i,i,t_huen(i+1));
    fprintf('\tk1=h*f(t_huen(%d),y1_huen(%d),y2_huen(%d))=%f\n',i,i,i,k1)
    fprintf('\tk1=h*g(t_huen(%d),y1_huen(%d),y2_huen(%d))=%f\n',i,i,i,l1)
  
    k2=h*f(t_huen(i+1),y1_huen(i)+k1,y2_huen(i)+l1);
    l2=h*g(t_huen(i+1),y1_huen(i)+k1,y2_huen(i)+l1);
    fprintf('\tk2=h*f(t_huen(%d+1),y1_huen(%d)+k1,y2_huen(%d)+l1)=%f\n',i,i,i,k2)
    fprintf('\tl2=h*g(t_huen(%d+1),y1_huen(%d)+k1,y2_huen(%d)+l1)=%f\n',i,i,i,l2)
  
    y1_huen(i+1)=y1_huen(i)+(1/2)*(k1+k2);
    y2_huen(i+1)=y2_huen(i)+(1/2)*(l1+l2);
    fprintf('\n\ty1_huen(%d+1)=y1_huen(%d)+(1/2)*(k1+k2)=%f\n',i,i,y1_huen(i+1))
    fprintf('\ty2_huen(%d+1)=y2_huen(%d)+(1/2)*(l1+l2)=%f\n',i,i,y2_huen(i+1))
end
fprintf('\n\t Using Modified Euler method at t=%f ,y(%f)=%f\n',t_huen(end),t_huen(end),y1_huen(end))
%%%%%%%%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%%%%%%%