0

I trying to make a Matlab code to plot a discrete solution of the heat equation using the implicit method.

The information I am given about the heat equation is the following:

d^2u/d^2x=du/dt

Initial conditions (t=0):

u=0 if x>0

Otherwise u=1 (when t=0)

The discrete implicit difference method can be written as follows:

(I+delta t*A)[v(m+1)]=v(m), where I is an identity matrix, delta t is the times space, m is the time-step number, v(m+1) is the v-value at the next time step. A is the matrix: A has the value 2 at the diagonal, while -1 both right below and right over this diagonal. All the other numbers are zero.

A=(1/dx^2)*
     [2 -1  0 . . . 0

    -1  2 -1 . . . .

     0 . . .  . .  .

     . . .  .   .  0

     . . . . -1  2 -1

     1 . . .  0 -1  2]

Without inserting the values, the symbols are written instead:

A=(1/dx^2)*
     [alpha(1) gamma(1)  0  . . .  . . . . . 0

      beta(2)  alpha(2) gamma(2) ... .. .   0

     0 . . .  . .  . . . . . . . . . . . . .

     . . .  .   . . . . . . . . . . . . .   0

     . . . .     beta(n-1)  alpha(n-1) gamma(n-1)

    0 . . .              0     beta(n)  alpha(n)]
L = 20.; % Length of the wire
T =0.1; % Final time
maxm = 2000; % Number of time steps
dt = T/maxm;
n = 70; % Number of space steps
dx = L/(n);
r = dt/(dx^2); % Stability parameter, r less or equal to 1/2



alpha(1)=2+(1/(r));
d(1)=alpha(1);
v(1)=0; % Boundary condition.
v(n)=1; %Boundary condition.
c(1)=v(0);

for k=2:n
    beta(k)=-1;
    gamma(k)=-1;
    alpha(k)=2+(1/r);
    m(k)=beta(k)/d(k-1);
    d(k)=alpha(k)-(m(k)*gamma(k-1));
    c(k)=v(k)-(m(k)*c(k-1));
    v(n)=c(n)/d(n);
    for k=n-1:1
        v(k)=(c(k)-gamma(k)*v(k+1))/d(k);
    end
end

The error message I get is the following:

Attempted to access v(2); index out of bounds because numel(v)=1.

Error in Implicit (line 37) c(k)=v(k)-(m(k)*c(k-1));

Does anyone know what I should write instead so that the error message disappear? Does the Matlab code look okay, or is there something I should change? David

David
  • 159
  • 8
  • You don't say in your code how `b` is defined, but the error message is pretty clear: `b` is a scalar (dimension 1) and you are attempting to access `b(2)` in your loop, hence the error. Up to you to fix it. – am304 Mar 09 '15 at 12:41
  • Thank you also for your reply. I do not know how I should define b. Do know how it should be defined so that it satifies the information at the top of the question? – David Mar 12 '15 at 09:48
  • You are trying to use `v(2)` when writing `c(k) = ...` but it's not defined until the following line `v(n) = c(n)/d(n)`, hence the error. `v` and `c` are mutually dependent, so you probably need to re-arrange your equations so that you can compute one, then the other. – am304 Mar 12 '15 at 09:54
  • Okay, thank you. I will try to re-arrange the code. – David Mar 12 '15 at 11:29

1 Answers1

0

I don't think your statement of initial conditions is correct.

Your equation has been normalized, obviously. Your body sounds like a semi-infinite solid. I'd write that initial condition this way: 

u(x, t) = 1 for t = 0

This says that the body is at a uniform temperature when the problem starts. For t > 0 you've got a boundary condition for temperature, flux, or convection at the surface that will remove energy from the body and give you a transient temperature distribution that will be of some interest.

duffymo
  • 305,152
  • 44
  • 369
  • 561
  • Thank you for your reply. I have made some changes in the code. This time I tried without using the initial condition, but only boundary conditions. The problem is that I do not know how I should define v(k) in order for the error message to disappear. Do you know how I can define v(k) from the information on the top of the text? – David Mar 11 '15 at 10:48
  • You use the same index in the inner and outer loops. That seems obviously bad to me. Since this is a matrix, I'd expect two indexes: one for rows, another for columns. Make them (i, j) or (m, n), or anyting else that makes sense to you. – duffymo Mar 11 '15 at 12:28
  • Okay, thank you again. The reason I use the same index in the inner an outer loops is because v(k+1) are to be found using backwards tracking. I only used what it said in my book. To me, it seems like the biggest problem is that I do not have v(2), v(3),..v(n-1). I only have v(1) and v(n). v(2), v(3),..,v(n-1) should have been used to find the v-values in the last for loop. Do you know what I can do in order to find v(2), v(3),..v(n-1)? I am only given the information at the top of the page. ( The v-values in the last for loop are v-values at the next time step). – David Mar 12 '15 at 08:08