38 lines
1.0 KiB
Matlab
38 lines
1.0 KiB
Matlab
function testRK4Programs
|
|
figure
|
|
hold on
|
|
t = 0:0.1:10;
|
|
y = 0.0625.*(t.^2+4).^2;
|
|
plot(t, y, '-k')
|
|
[tode4, yode4] = testODE4(t);
|
|
plot(tode4, yode4, '--b')
|
|
[trk4, yrk4] = testRK4(t);
|
|
plot(trk4, yrk4, ':r')
|
|
legend('Exact', 'ODE4', 'RK4')
|
|
hold off
|
|
fprintf('Time\tExactVal\tODE4Val\tODE4Error\tRK4Val\tRK4Error\n')
|
|
for k = 1:10:length(t)
|
|
fprintf('%.f\t\t%7.3f\t\t%7.3f\t%7.3g\t%7.3f\t%7.3g\n', t(k), y(k), ...
|
|
yode4(k), abs(y(k)-yode4(k)), yrk4(k), abs(y(k)-yrk4(k)))
|
|
end
|
|
end
|
|
|
|
function [t, y] = testODE4(t)
|
|
y0 = 1;
|
|
y = ode4(@(tVal,yVal)tVal*sqrt(yVal), t, y0);
|
|
end
|
|
|
|
function [t, y] = testRK4(t)
|
|
dydt = @(tVal,yVal)tVal*sqrt(yVal);
|
|
y = zeros(size(t));
|
|
y(1) = 1;
|
|
for k = 1:length(t)-1
|
|
dt = t(k+1)-t(k);
|
|
dy1 = dt*dydt(t(k), y(k));
|
|
dy2 = dt*dydt(t(k)+0.5*dt, y(k)+0.5*dy1);
|
|
dy3 = dt*dydt(t(k)+0.5*dt, y(k)+0.5*dy2);
|
|
dy4 = dt*dydt(t(k)+dt, y(k)+dy3);
|
|
y(k+1) = y(k)+(dy1+2*dy2+2*dy3+dy4)/6;
|
|
end
|
|
end
|