Midpoint numerical integration without a built in function

268 views (last 30 days)
Aggie
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.
thanks in advance.
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);
  1 Comment
Roger Stafford
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.

Sign in to comment.

Answers (1)

Geoff Hayes
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)));
So your code becomes
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?

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!