Need help with Runge Kutta 4th order code

Question

Peter Phung 2018-10-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/424367-need-help-with-runge-kutta-4th-order-code

9780470501979.pdf

I am trying to plot example 5.3 in my heat transfer textbook using 4th order Runge Kutta and on the graph in the text book, it says that tc happens at 124 seconds. Unfortunately for my code, it shows that tc happens at 123 seconds. If I shift the time by one place, it will make tc = 124 seconds, but then the initial temperature (25 degrees celsius) will be shifted to 1 second instead of 0 seconds which is not what the graph in the textbook depicts. tc is defined as the time that it takes for the temperature to reach 150 seconds. I have attached the PDF of the question and my code is provided below.

clear all; close all; clc
tstart=0;
tend=150;
h = 1;   %%change to small step size (0.1 or 0.01) then the result of numerical method will be better
t = tstart:h:tend;
w(1) = 25+273;
t(1) = 0;
T(1) = 25+273;
%fprintf('Step 0: t = %12.8f, w = %12.8f\n', t, w);
for i=2:size(t,2)
k1 = h*f(t(i-1),w(i-1));
k2 = h*f(t(i-1)+h/2, w(i-1)+k1/2);
k3 = h*f(t(i-1)+h/2, w(i-1)+k2/2);
k4 = h*f(t(i-1)+h, w(i-1)+k3);
w(i) = w(i-1) + (k1+2*k2+2*k3+k4)/6;
end
figure (1)
neww = w-273 %subtract 723 kelvins from w matrix to convert kelvins back into degrees
plot(t,neww)
title('Plot of temperature vs. time for example 5.3 (first 150 seconds)')
xlabel('Time (seconds)')
ylabel('Temperature (degrees)')
%%%%%%%%%%%%%%%%%%
function deltaT = f(t,T)
convection = 40
tsurr = 175+273
tinf = 175+273
emissivity = 0.8
const = 5.67*10^-8
roe = 2770
ccc = 875
length = (3*10^-3)/2
deltaT = -(convection*(T-tinf)+(emissivity*const*((T^4)-(tsurr^4))))/(roe*ccc*length);
end