Re: Modelling questions

[c.]

On 23 Jun 2011, at 15:55, И Страшко wrote:

> Sorry about my last post.  Question 1 was really whether this way of 
> iterating was applicable.  I forgot to ask that, and pushed "send" before 
> realizing I didn't ask it.  
> 
> Yes, I used fsolve.  Let me ask my question more simply.
> 
> 
> I took my pde, for which dt (really deltat) and T_e are constants,
> 
> k * del2 T = dT/dt - T_e, and rewrote it in a form,
> 
> 0 = k * del2(T) - dT/dt + T_e
> 
> dT/dt is approximated using 
> 
> dT/dt = (T(i) - T(i-1)) / dt, which is a kind of forward derivative.
> 
> so
> 
> 0 = k * del2(T(i)) - (T(i) - T(i-1)) / dt + T_e
> 
> 
> 
> Next, I prepared a function, called heat.  The variable x here is T(i) above,
> the variable Tlast is T(i-1), and Te is T_e.  
> 
> function y=heat(x)
>   global K dx dy dz dt Te Tlast
>   [xl,yl,zl]=size(x);
>   z = K*del2(x,dx,dy,dz);
>   z = z + Te*ones(xl,yl,zl);
>   z = z-(x-Tlast)/dt;
>   y = z;
> endfunction
> 
> I then use this in fsolve.  
> 
> Then I advance from one time to the next.  It seems to work, it looks 
> physical, and it converges well.  
> 
> The question is... Is this mathematically correct?  (BTW, why should one 
> write a second order routine when Octave supplies del2?)
> 
> 
> (My second question seems to be an adaptation of the FEM, so I will pursue 
> that analysis as described in Wikipedia.)
> 
> Thanks for the help.
> 
> -John Stasko

If I understand correctly you wish to solve the nonstationary heat eqaution 
using centered finite differences for space discredization and an implict time 
stepping algorithm.

You don't need to do the time advancement yourself as Octave can take care of 
that for you with lsode or daspk, actually
with the package BIM it can take care of the space discretization as well, but 
that is trivial to do manually if you don't 
want to install an additional package.

pkg load bim

N  = 11;                  % number of space discretization points
M  = 11;                  % number of time discretization points

x  = linspace (0, 1, N).';  % space discretization grid
t  = linspace (0, 1, M).';  % space discretization grid

k  = 1;                   % diffusion coefficient
Te = 1;                   % source term NB with these coefficients 
                          % the solution will never reach a stationary value

u0 = zeros (N, 1);        % initial condition

A  = bim1a_laplacian (x, k * ones (N-1, 1), ones (N, 1)); % discretization of - 
k Laplacian 
M  = bim1a_reaction  (x, ones (N-1, 1), ones (N, 1));     % mass matrix
b  = bim1a_rhs       (x, ones (N-1, 1), Te *ones (N, 1)); % source term

% your time discretized PDE is expressed now as 
%     res (u, udot, t) =  M * udot + A * u - b = 0

res = @(u, udot, t) M * udot + A * u - b;

% DASPK will solve this equation in space-time for you
u = daspk (res, u0, u0, t);

surf (x, t, u)
xlabel ('x')
xlabel ('t')

HTH,
c.