Question

Solve this differential equation using Matlab

yy' + xy² =x , with y(0)=5 for x=0 to 2.5 with a step size 0.25

(a) Analytical

(b) Euler

(c) Heun

d) 4th order R-K method

Display all results on the same graph

Answer 1

MATLAB Code:

close all
clear
clc

dydx = @(x,y) (x - x*y^2)/y; % Given ODE, y'
h = 0.25; % Step size
xi = 0; xf = 2.5;
x = xi:h:xf;
n = (xf - xi)/h; % Number of nodes

% Part (a)
% Analytical solution
syms y_exact(x_exact)
eqn = y_exact*diff(y_exact, x_exact) + x_exact*y_exact^2 == x_exact; % ODE
cond = y_exact(0) == 5; % Initial condition
y_exact = dsolve(eqn, cond); % Solve for exact solution
x_vals = linspace(xi, xf); % Bunch of new samples for a smoother plot
plot(x_vals, subs(y_exact, x_vals)), hold on

% Part (b)
% Euler's Method
y_euler(1) = 5; % Initial condition
for i=1:n
m = dydx(x(i), y_euler(i));
y_euler(i+1) = y_euler(i) + h*m; % Euler's Update
end
plot(x, y_euler, 'o-')

% Part (c)
y_heun(1) = 5; % Initial condition
for i = 1:n
m1 = dydx(x(i), y_heun(i));
m2 = dydx(x(i) + h, y_heun(i) + h*m1);
y_heun(i+1) = y_heun(i) + 0.5*h*(m1 + m2); % Heun's Update
end
plot(x, y_heun, '*-')

% Part (d)
% RK4 Method
y_rk4(1) = 5; % Initial condition
for i = 1:n
k1 = dydx(x(i), y_rk4(i));
k2 = dydx(x(i) + 0.5*h, y_rk4(i) + 0.5*h*k1);
k3 = dydx(x(i) + 0.5*h, y_rk4(i) + 0.5*h*k2);
k4 = dydx(x(i) + h, y_rk4(i) + k3*h);
y_rk4(i + 1) = y_rk4(i) + (1/6)*(k1 + 2*k2 + 2*k3 + k4)*h; % RK4 Update
end
plot(x, y_rk4, 's-')

xlabel('x'), ylabel('y(x)'),
title('Solution of ODE')
legend('Exact Solution', 'Euler', 'Heun', 'RK4')

Plot: