[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Modelling questions

From: Thomas Shores
Subject: Re: Modelling questions
Date: Sun, 26 Jun 2011 02:59:46 -0500
User-agent: Mozilla/5.0 (X11; U; Linux x86_64; en-US; rv: Gecko/20110428 Fedora/3.1.10-1.fc14 Thunderbird/3.1.10

On 06/23/2011 03:31 PM, c. wrote:
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.


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
   z = K*del2(x,dx,dy,dz);
   z = z + Te*ones(xl,yl,zl);
   z = z-(x-Tlast)/dt;
   y = z;

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')

Help-octave mailing list
No, I do not think that your approach -- at least as you have described it in two posts -- is mathematically correct. As Carlo suggests, you seem to be approximating the PDE by implicit time stepping (backwards Euler) and centered central differences for second derivatives. If correctly formulated, this method is well known to be unconditionally stable and first order convergent in time step, second order in space step. However, you seem to largely ignore boundary conditions, except that you force them at the end of each iteration. Boundary conditions will affect the central second differences next to the boundary, and you do not account for it. In fact, you seem to assume that the PDE holds *at* the boundary, which need not be correct. Failure to account for these facts will likely make your method at best conditionally stable and it will certainly not realize the accuracy of the same method developed with correct attention to the boundaries.

A common technique in numerical analysis is to try a method on a test problem with a known solution so that you can see the error behavior. For example, you might start with the solution
T(x,t) = exp(-t)*sin(x) +x*(pi-x), 0<=x<=pi, 0<=t,
which satisfies the diffusion equation with source term Te = 2 and diffusion coefficient K = 1 on the space-time intervals indicated. Try your method with this problem and see how refining step size in time and space works as regards error reduction. If you have further comments or questions relating to your code, I suggest that you explicitly list the code in question.

Thomas Shores

reply via email to

[Prev in Thread] Current Thread [Next in Thread]