Hey everyone,(adsbygoogle = window.adsbygoogle || []).push({});

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:

http://img165.imageshack.us/img165/3875/2ndorderodetm7.jpg [Broken]

Which I've rearranged to:

d^{2}y_{2}/dx^{2}= phi1 * y_{1}+ phi2 * (y_{1})^{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...

%ThieleBVP.m - Example 4

clear

clc

global phi

global phi2

global constant

%Maximum number of secant iterations

maxsec = 20;

%Absolute error tolerance

tol = 1e-15;

%Arrays to store the guessed boundary condition, c(0), and the residual, R

c0 = zeros(maxsec, 1);

R = zeros(maxsec, 1);

%Set the phi is for T, phi2 is for T^4, constant is the rest of the

%collected terms

%Setting Constants...

kcon = 72

h = 2000

eps = 0.1

sig = 5.67 * 10^-8

Icur = 2

rowe = 32*10^-8

Tinf = 300

D = 80*10^-6

L = 4*10^-3

%Calculating Stuff...

RHS = -(Icur^2*rowe)/(kcon*(pi/4*D^2)^2)

LHS1 = -4*h/(kcon*D)

LHS2 = -4*eps*sig/(kcon*D)

%Finding the final constants...

phi = -LHS1

phi2 = -LHS2

constant = RHS + LHS1*Tinf + LHS2*Tinf

%Target boundary condition to shoot for, c(1)

c1_true = 300;

%First guess of the unknown boundary condition, c0 = c(0)

c0(1) = 280;

%Shooting loop

%rspan = [0 1]; % integration interval

rspan = [-L/2 L/2];

% first shot

i = 1;

yini = [c0(i) 0]; % initial conditions

[r,y] = ode45('ODEs',rspan,yini);

R(i) = y(length(r),1) - c1_true;

%second shot

i = 2;

c0(i) = 1.05*c0(i-1);

yini = [c0(i) 0]; % initial conditions

[r,y] = ode45('ODEs',rspan,yini);

R(i) = y(length(r),1) - c1_true;

while abs(R(i))>tol

i = i+1;

c0(i) = c0(i-1)-R(i-1)*(c0(i-1)-c0(i-2))/(R(i-1)-R(i-2));

yini = [c0(i) 0]; % initial conditions

[r,y] = ode45('ODEs',rspan,yini);

R(i) = y(length(r),1) - c1_true;

if i == maxsec

display('Maximum no. of iterations reached!');

break;

end

end

plot(r,y(:,1),'ro');This just results in MATLAB spitting out jibberish. function dydr = ODEs(r,y)

global phi

global phi2

global constant

dydr = zeros(2,1);

dydr(1) = y(2);

dydr(2) = phi*y(1)+phi2*(y(1))^4-constant;

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:

http://img379.imageshack.us/img379/6339/80553833xt7.jpg [Broken]

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...

**Physics Forums - The Fusion of Science and Community**

Dismiss Notice

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# 2nd Order Non-Linear ODE in MATLAB Issues

Loading...

Similar Threads - Order Linear MATLAB | Date |
---|---|

Matlab Simulink: designing a 2nd order sliding controller | Jul 5, 2017 |

LaTeX Changing the order of references using bibtex to sequential | Mar 8, 2016 |

Matlab help - 2nd order ODE | Mar 12, 2015 |

[Matlab] Order of precedence in computations | Jul 30, 2014 |

**Physics Forums - The Fusion of Science and Community**