Path: news.mathworks.com!not-for-mail
From: "Demoslav " <dusanori@gmail.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Finite difference method (implicit) for 1-D heat equation
Date: Wed, 14 Dec 2011 00:41:08 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 43
Message-ID: <jc8rb4$m87$1@newscl01ah.mathworks.com>
References: <jc8qf1$jqc$1@newscl01ah.mathworks.com>
Reply-To: "Demoslav " <dusanori@gmail.com>
NNTP-Posting-Host: www-05-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1323823268 22791 172.30.248.37 (14 Dec 2011 00:41:08 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Wed, 14 Dec 2011 00:41:08 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 3227944
Xref: news.mathworks.com comp.soft-sys.matlab:752280

"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)...
> Please help I need this until thursday. Thanks for your help.
> 
> 
> 
> function [U]=HE1D(a2, x0, xn, t0, tm, h, k, Tzac, Xz, Xk)
> % equation: a2*u_t=u_xx , (x0, xn)x(t0, tm) 
> % Tzac - initial condition, Xz,Xk - boundary conditions
> % exact solution u = exp(t-x)
> %[U]=HE1D(1, 0, 1, 0, 1, 0.01, 0.01, 'exp(-x)' , 'exp(t)', 'exp(t-1)')
> 
> a=xn-x0; b=tm-t0;           
> n=floor(a/h); m=floor(b/k);         
> 
> % initial and boundary conditions
> zac=inline(Tzac);
> okr1=inline(Xz);
> okr2=inline(Xk);
> t=[t0+k:k:tm];
> x=[x0:h:xn];
> 
> % solution matrix initialization
> U=ones(m+1, n+1); 
> U(1, :)=zac(x);
> U(2:m+1, 1)=okr1(t); 
> U(2:m+1, n+1)=okr2(t);
> 
> % the 3-daig matrix
> R=k/(a2*h^2);
> D=1+2*R;
> P=ones(n-1,1);
> M=diag(D*ones(n-1,1)) - diag(R*ones(n-2,1),1) - diag(R*ones(n-2,1),-1);
> 
> for i=2:m+1
>     P(1)=U(i-1,2)+R*U(i,1);
>     P(end)=U(i-1,n)+R*U(i,n+1);
>     P(2:end-1)=U(i-1,3:n-1);
>     u=M\P;
>     U(i,2:end-1)=u;
> end

please close/delete/ignore the my first thread, thanks.