Question

Solve using matlab code

The initial value problem dydx−y= 2 cosx, y(0) =−2

has the exact solution y(x) =−e^x −√2 cos (x+π4).

Use the Euler method to solve the initial value problem for 0≤x≤2 using n=10,50,100,200 and plot solutions in one graph.

Repeat #1 using the Runge-Kutta method and plot solutions in one graph with the exact solution

Answer 1

MATLAB Script:

close all
clear
clc

% Given
f = @(x,y) 2*cos(x) + y; % Given ODE
x0 = 0; xf = 2; % Intervals of x
y0 = -2; % Initial condition
y_exact = @(x) -exp(x) - sqrt(2)*cos(x + pi/4); % Exact solution

% Euler Method Solutions
n1 = 10;
h1 = (xf - x0)/n1; % Step Size 1
x1 = x0:h1:xf;
y1 = my_euler(x0, y0, xf, h1, f);

n2 = 50;
h2 = (xf - x0)/n2; % Step Size 2
x2 = x0:h2:xf;
y2 = my_euler(x0, y0, xf, h2, f);

n3 = 100;
h3 = (xf - x0)/n3; % Step Size 3
x3 = x0:h3:xf;
y3 = my_euler(x0, y0, xf, h3, f);

n4 = 200;
h4 = (xf - x0)/n4; % Step Size 4
x4 = x0:h4:xf;
y4 = my_euler(x0, y0, xf, h4, f);

% Plotting
figure, plot(x4, y_exact(x4)), hold on % Plot exact solution
plot(x1, y1, x2, y2, x3, y3, x4, y4), hold off % Plot Euler Method solutions
xlabel('x'), ylabel('y(x)'), title('Euler''s Method')
legend('Exact Solution', 'Euler''s Method (n = 10)', ...
'Euler''s Method (n = 50)', 'Euler''s Method (n = 100)', ...
'Euler''s Method (n = 200)', 'Location', 'southwest')

% RK4 Method Solutions
n1 = 10;
h1 = (xf - x0)/n1; % Step Size 1
x1 = x0:h1:xf;
y1 = my_rk4(x0, y0, xf, h1, f);

n2 = 50;
h2 = (xf - x0)/n2; % Step Size 2
x2 = x0:h2:xf;
y2 = my_rk4(x0, y0, xf, h2, f);

n3 = 100;
h3 = (xf - x0)/n3; % Step Size 3
x3 = x0:h3:xf;
y3 = my_rk4(x0, y0, xf, h3, f);

n4 = 200;
h4 = (xf - x0)/n4; % Step Size 4
x4 = x0:h4:xf;
y4 = my_rk4(x0, y0, xf, h4, f);

% Plotting
figure, plot(x4, y_exact(x4)), hold on % Plot exact solution
plot(x1, y1, x2, y2, x3, y3, x4, y4), hold off % Plot RK4 Method solutions
xlabel('x'), ylabel('y(x)'), title('RK4 Method')
legend('Exact Solution', 'RK4 Method (n = 10)', ...
'RK4 Method (n = 50)', 'RK4 Method (n = 100)', ...
'RK4 Method (n = 200)', 'Location', 'southwest')

function y = my_euler(t0, y0, tf, h, f)
y(1) = y0;
t = t0:h:tf;
for i = 1:length(t)-1
y(i+1) = y(i) + h*f(t(i), y(i)); % Euler Update
end
end

function y = my_rk4(t0, y0, tf, h, f)
y(1) = y0;
t = t0:h:tf;
for i = 1:length(t)-1
k1 = f(t(i), y(i));
k2 = f(t(i) + 0.5*h, y(i) + 0.5*h*k1);
k3 = f(t(i) + 0.5*h, y(i) + 0.5*h*k2);
k4 = f(t(i) + h, y(i) + k3*h);
y(i + 1) = y(i) + (1/6)*(k1 + 2*k2 + 2*k3 + k4)*h; % RK4 Update
end
end

Plots: