Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Rugga Kutta won't converge for 2nd order ODE

Subject: Rugga Kutta won't converge for 2nd order ODE

From: Windell

Date: 8 Dec, 2009 11:32:03

Message: 1 of 2

Hi,

I've written a program in MATLAB that uses the Rugga Kutta ode45 solver to interatively solve for the temperature distrubution in a cylinder undergoing 1-dimensional 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.715287e-002. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.110223e-016) at time t.

Can anyone give suggestions? I know that my differential equations are correct. I've also reduced my tolerances to 1*e-3. I also tried the other ODE solvers such as ode23,ode113, etc. They all fail.

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',1e-3,'AbsTol',1e-3)
%call the solver.
[x,p]=ode45('temp_distribution2',[0 0.1],p0);

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;
%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)^4-T_i^4)+h*(p(1)-T_i));

Subject: Rugga Kutta won't converge for 2nd order ODE

From: Steven Lord

Date: 8 Dec, 2009 14:55:06

Message: 2 of 2


"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 Runge-Kutta.

http://en.wikipedia.org/wiki/Runge-Kutta

> to interatively solve for the temperature distrubution in a cylinder
> undergoing 1-dimensional 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.715287e-002. Unable to meet integration tolerances without
> reducing the step size below the smallest value allowed (1.110223e-016) at
> time t.
>
> Can anyone give suggestions? I know that my differential equations are
> correct. I've also reduced my tolerances to 1*e-3. 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.715287e-002, 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',1e-3,'AbsTol',1e-3)

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)^4-T_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.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us