# Midpoint numerical integration without a built in function

268 views (last 30 days)
Aggie on 20 Oct 2014
Commented: Roger Stafford on 20 Oct 2014
I need some help building a matlab script to solve dy/dt = y*t^3-1.5*y using the midpoint method. I have solved this using Euler's and the below code
a = 0;
b = 2;
h=.5;
T = 0:h:2;
y = zeros(1,((b-a)/h+1));
y(1) = 1;
phi1 = y.*T.^3-1.5.*y;
for i = 2:length(T)
phi1 = y(i-1).*T.^3-1.5.*y(i-1);
y(i) = y(i-1) + phi1(i-1).*h;
end
plot(T,y,':bo')
But solving cannot figure out the midpt method as I know the +1/2 intervals are tough on MATLAB.
Below is what I have for midpoint and I know it is very wrong.
a = 0;
b = 2;
h= .5;
h2 = h/2;
t = 0:h2:2;
y = zeros(1,9);
y(2) = 1;
phi3 = y.*t.^3-1.5.*y;
for i = 3:(length(t))
Y2 = y(i-2)+ (y(i-2).*(t).^3-1.5.*y(i-2))*(h/2);
phi3 = Y2.*t.^3-1.5.*Y2;
y(i) = y(i-2) + phi3(i-1).*h;
end
plot(t,y);
Roger Stafford on 20 Oct 2014
It pains me to see you use numerical integration for this differential equation when its solution can be expressed so easily as
y = K*exp(t^4/4-1.5*t)
where K depends on initial conditions.

Geoff Hayes on 20 Oct 2014
Aggie - the midpoint method should be very similar to your Euler implementation, with just a couple of minor changes (for example the step size). The above code gets a little confusing with the re-calculation of phi3 at every iteration (as it is vector of which you only use one element from), so you may want to try an alternative approach where you define a function handle to your equation and evaluate that instead
F = @(t,y)y*t^3-1.5*y;
Then, the difference relation for the midpoint algorithm can be written as
y(n) = y(n-1) + h*F(t(n-1) + h/2, y(n-1) + (h/2)*F(t(n-1),y(n-1)));
a = 0;
b = 2;
h = 0.5;
F = @(t,y)y*t^3-1.5*y;
t = a:h/2:b;
y = zeros(size(t));
% set the initial condition
y(1) = 1;
for n=2:length(t)
y(n) = y(n-1) + h*F(t(n-1) + h/2, y(n-1) + (h/2)*F(t(n-1),y(n-1)));
end
plot(t,y,':ro');
Note that the body of the for loop has been reduced to one line, and is quite different from yours. Is there a reason that you started iterating at i equal to 3, and referred to i-2 in your equation?