Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Finite difference method (implicit) for 1-D heat equation
Date: Wed, 14 Dec 2011 05:17:07 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 30
Message-ID: <jc9bgj$8b5$1@newscl01ah.mathworks.com>
References: <jc8qf1$jqc$1@newscl01ah.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: www-01-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1323839827 8549 172.30.248.46 (14 Dec 2011 05:17:07 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Wed, 14 Dec 2011 05:17:07 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: news.mathworks.com comp.soft-sys.matlab:752290

"Demoslav " <dusanori@gmail.com> wrote in message <jc8qf1$jqc$1@newscl01ah.mathworks.com>...
> I have a problem with my code. Im trying to solve the 1-D heat equation via implicit finite difference method. Ive been looking at my code for hours but I cant find no mistake. Problem is that I get huge errors, I mean when I substract the exact solution from the numerical it is much bigger then o(h*h + k)...
> .......
- - - - - - - - - -
  As a check on your procedure I would suggest that you sample your matrix U at various points (i,j) in its interior and its edges to verify that the quantities

 (U(i-1,j-1)-2*U(i-1,j)+U(i-1,j+1))*h^2

and

 a2*(U(i,j)-U(i-1,j))*k

are indeed equal (up to round off errors.)  In fact you could easily write a short routine that would find the maximum absolute deviation between them for the entire matrix U.  These are the finite difference equations that you are presumably finding the solution to as an approximation to the partial differential equation a2*dU/dt = d^2U/dx^2.

  I would also try using smaller values for h and k to increase the accuracy of the above difference equations as compared with the ideal partial differential equations to see if that brings you significantly closer to the precise solution.

  Finally, I would also suggest that you modify the two lines

 a=xn-x0; b=tm-t0;           
 n=floor(a/h); m=floor(b/k);         

to

 a=xn-x0; b=tm-t0;           
 n=round(a/h); m=round(b/k); 
 h = a*n; k = b*m;

to avoid possible discrepancy between h and k and the corresponding integers n and m.        

Roger Stafford