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.
_______________________________________________
Help-octave mailing list
address@hidden
https://mailman.cae.wisc.edu/listinfo/help-octave