- #1
TurboRegal
- 2
- 0
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:
http://img165.imageshack.us/img165/3875/2ndorderodetm7.jpg [Broken]
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:
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...
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:
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...
%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');
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;
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:
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...
Last edited by a moderator: