Midpoint numerical integration without a built in function

59 次查看（过去 30 天）
Aggie
Aggie 2014-10-20
评论： Konner Brickey 2022-10-28
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 个评论
Roger Stafford
Roger Stafford 2014-10-20
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Geoff Hayes
Geoff Hayes 2014-10-20
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?

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by