Question

10.16: Write a user-defined MATLAB function that solves a first-order ODE by applying the midpoint method (use the form of second-order Runge-Kutta method, Eqs(10.65),(10.66)). For function name and arguments use [x,y]=odeMIDPOINT(ODE,a,b,h,yINI). The input argument ODE is a name for the function that calculates dy/dx. It is a dummy name for the function that is imported into odeMIDPOINT. The arguments a and b define the domain of the solution, h is step size; yINI is initial value. The output arguments, x and y, are vectors with the x and y coordinates of the solution. Use the function odeMIDPOINT to solve the ODE in Problem 10.2. Write a MATLAB program in a script file that solves the ODE twice, once by using h=0.8 and once using h=0.1. The program should also plot exact solution (given in prob. 10.2 and two numerical solutions (all in same figure).10.2: first-order ODE: dy/dx=x-xy/2 from x=1 to x=3.4 with y(1)=1. Use h=0.8a) solve with Euler's explicit method.b) solve with modified Euler methodc) solve with classical fourth-order Runge-Kutta method. The analytical solution is y=2-e^((1-x^2)/4). In each part, calculate the error between the true solution and numerical solution at the points where the numerical solution is determined.

Answer 1

MATLAB Script:

close all
clear
clc

ODE = @(x,y) x - x*y/2;
a = 1; b = 3.4; yINI = 1;

h = 0.8;
[x1, y1] = odeMIDPOINT(ODE, a, b, h, yINI);

h = 0.1;
[x2, y2] = odeMIDPOINT(ODE, a, b, h, yINI);

xx = a:0.01:b;
plot(xx, 2 - exp((1 - xx.^2)/4), 'linewidth', 2), hold on
plot(x1, y1, 'o-', x2, y2, 'o-'), xlabel('x'), ylabel('y')
legend('Exact Solution', 'h = 0.8', 'h = 0.1', 'location', 'northwest')

function [x, y] = odeMIDPOINT(ODE, a, b, h, yINI)
x = a:h:b;
y(1) = yINI;
for i = 1:length(x)-1
k1 = h*ODE(x(i), y(i));
k2 = h*ODE(x(i) + h/2, y(i) + k1/2);
y(i+1) = y(i) + k2;
end
end

Plot: