Potential in a non uniform resistor.

My apologies. Google Groups seems to have wrapped my posting badly. Here is the same post with better formatting:

I need to calculate potential in a nonuniform conductor. The complete problem is 2 dimensional, but even 1D is causing me trouble.

In the 1D problem, you have a 1D conductor of length 3l. The first l has a conductivity of k and the remaining 2l of the conductor has a conductivity of 2k. The two sections have the same resistance. If a voltage is applied across this conductor, the both of the sections will drop half the voltage, and the voltage change in each section will be linear.

Here are the the equations:

Let J be current density, E be the electric field and V be the potential.

Since there is no free charge:

div J = 0.

The current is related to the field via the isotropic but not constant conductivity k:

J = k E.

And the field is the gradient of the potential V:

E = grad V.

Putting this together gives:

div (k grad V) = 0.

In 1D (dimension x) this is:

dk / dx * dv/dx + k * d^2 v / dx^2 = 0

d^2 v / dx^2 = - dk / dx * dv/dx

Taking a finite difference approximation gives:

k(i) * ( v(i-1) -2 v(i) + v(i+1) ) / D^2 = -(k(i+1) - k(i-1))/2D * (v(i+1) - v(i-1))/2D

where D is the grid spacing, and i is the index which is the descrete position.

This comes out as:

v(i) = = (v(i+1) * ( k(i) + (k(i+1) - k(i-1))/4) + v(i-1)*(k(i) - (k(i+1)-k(i-1))/4))/ (2 * k(i))

which can be iterated. This is the relaxation solution for the differential equation above.

If I solve the above equation, for the problem above, I seem to get the right answer (the middle voltage gets very close to half the total voltage, butit takes a long time to converge). However, if I make the conductivitites significantly different, then the equations become unstable and I can't get an answer.

Can anyone help?

-Ed

I have a simple solver in matlab. Here is the code:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% xs = [1 : 90 ]; %Positions vs = 0 * xs; %Voltages

rs = 0 * xs; %Conductivities

rs(1:30) = 5; rs(31:90) = 10;

while 1

vs(1) = 1; %Boundary conditions: 1 volt at the left edge vs(90) = 0; %Boundary conditions: 0 volts at the right edge edge

v2 = vs;

for i=2:89 v2(i) = (vs(i+1) * ( rs(i) + (rs(i+1) - rs(i-1))/4) + ... vs(i-1)*(rs(i) - (rs(i+1)-rs(i-1))/4))/ (2 * rs(i)); end

hold off plot(xs, vs, 'g'); hold on plot(xs, v2, 'r');

vs = v2;

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%