"Windell " <windelljones@gmail.com> wrote in message
news:hfldfj$lba$1@fred.mathworks.com...
> Hi,
>
> I've written a program in MATLAB that uses the Rugga Kutta ode45 solver
FYI, the name of the solver is RungeKutta.
http://en.wikipedia.org/wiki/RungeKutta
> to interatively solve for the temperature distrubution in a cylinder
> undergoing 1dimensional heat transfer with radiation, convection, and
> conduction. However, when I run the program, the Rugga Kutta fails giving
> me the error:
>
> Failure at t=4.715287e002. Unable to meet integration tolerances without
> reducing the step size below the smallest value allowed (1.110223e016) at
> time t.
>
> Can anyone give suggestions? I know that my differential equations are
> correct. I've also reduced my tolerances to 1*e3. I also tried the other
> ODE solvers such as ode23,ode113, etc. They all fail.
Run your code and then plot x versus the columns of p. You should note that
p makes a SHARP climb at about x = 4.715287e002, but even before then the
function value grows very quickly.
> Any suggestions?
> Thank you for your time,
> Windell
>
> Here is my code:
> %define constants.
> V=15.0;I=0.55;Tb=87.9+273.15;
> %calculate the input power at the base of the rod.
> P=V*I;
> %define constants.
> k=121; D=0.01; A=(3.14159/4)*D^2;
> %define initial conditions.
> p0=[Tb P/(k*A)];
> options = odeset('RelTol',1e3,'AbsTol',1e3)
This value for AbsTol is unreasonably low. Look at the values in your
initial condition vector, and look at the values in the p matrix
(particularly towards the end of the run.)
> %call the solver.
> [x,p]=ode45('temp_distribution2',[0 0.1],p0);
Note that ODESET does ***NOT*** make "global changes" to the ODE solver
state  in order for the tolerance changes you made in your call to ODESET
to apply here, you MUST pass the options structure into the ODE45 call.
Even making that modification didn't help matters.
> function dpdx=temp_distribution2(x,p)
> %Define constants.
> k=121.0; D=0.01; E=1.0; sigma=5.670*10^(8);A=(3.14159/4)*D^2;
> %state the temperature of the room (large enclosure approximation)
> T_i=21.9+273.15;
I think you're going to HAVE to rescale this problem. Note that in your
expression for dpdx(2) you're raising this quantity to the _fourth power_ 
that means that expression includes a value on the order of 1e9!
> %define system of diffeqs.
> dpdx=zeros(size(p));
> dpdx(1)=p(2);
> h=(0.36*k/D);
> Per=3.14159*D;
> dpdx(2)=(Per/k/A)*(E*sigma*(p(1)^4T_i^4)+h*(p(1)T_i));
So you're multiplying a very small number (relatively speaking), sigma,
times the difference of two very large numbers. This screams "catastrophic
cancellation" to me.
http://en.wikipedia.org/wiki/Loss_of_significance
I think you should go back and change the units in your problem to avoid
working with such large and small numbers. So instead of working with
temperatures in Kelvin (as I suspect you're doing, from the fact that your
constant T_i involves the quantity 273.15) work in Celcius (T_i2 = T_i 
273.15) and make the necessary adjustments to your equations.

Steve Lord
slord@mathworks.com
comp.softsys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
