2nd Order Non-Linear ODE in MATLAB Issues

1. Jul 10, 2008

TurboRegal

Hey everyone,
Having some trouble here using the solver we were supplied and modifying it to fit our problem...

I have a wire with a current flowing through it. I'm trying to find the temperature distribution wrt. position in the wire.

BV's are:
T(x=L/2) = 300K
dT/dx (x=0) = 0
(Apparently this is suppose to result in a symmetrical distribution around 0 from -L/2 -> L/2)

The DE is as follows:

Which I've rearranged to:
d2y2/dx2 = phi1 * y1 + phi2 * (y1)4 - constants

Where y1 = T, y2 = dT/dx, and phi1, phi2, constants are just the collected constant terms...

Here is the code I'm using, in it c->T, r->x, sorry the terms were different and I haven't had time to change them yet...
This just results in MATLAB spitting out jibberish.

At the moment, I need to figure out the following:
1) Does ode45 work for nonlinear ODE's?
2) How do I set the dT/dx (x=0) = 0 boundary value? I know this isn't working because if I set the phi2 power above to ^2 instead of ^4, I get:

3) Why isn't MATLAB solving this properly?
4) Is the "solver" part of the 1st code even going to work for this type of problem. It was used at first for a 2nd order linear ODE, I'm assuming this is what is causing jibberish?

Thanks for the help...

Last edited: Jul 10, 2008
2. Jul 10, 2008

TurboRegal

Fixing a few small errors in my stupid code...

%Finding the final constants...
phi = -LHS1
phi2 = -LHS2
constant = RHS + LHS1*Tinf + LHS2*Tinf^4

I think that's all I changed...

Nets me:

And gives me the error:
Warning: Failure at t=-1.512856e-003. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.469447e-018) at time t.
> In ode45 at 355
In ThieleBVP at 70

When I try just going from 0->L/2 and then making is symmetric myself it nets me the following:

This nets me (NO ERRORS):

But it should be going up in the middle no?

3. Jul 11, 2008

Jitse Niesen

1. Yes

2. If you set TSPAN (the second argument of ode45) to [-L/2 L/2] then Matlab applies the initial conditions at -L/2 and integrates from x = -L/2 to x = L/2. Since you want to apply the initial conditions at x = 0, you have to set TSPAN to [0 L/2] (as you do below). Then you can compute the other half of the solution by symmetry, or by doing another integration but with TSPAN set to [0 -L/2].

What do you mean with it should be going up in the middle? The solution looks fine to me, except that you get a temperature of -600K at x=0 which is physically impossible. Perhaps make sure that your constants are correct?