Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Finite difference method (implicit) for 1-D heat equation

Subject: Finite difference method (implicit) for 1-D heat equation

From: Demoslav

Date: 14 Dec, 2011 00:26:09

Message: 1 of 5

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

Subject: Finite difference method (implicit) for 1-D heat equation

From: Demoslav

Date: 14 Dec, 2011 00:41:08

Message: 2 of 5

"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.

Subject: Finite difference method (implicit) for 1-D heat equation

From: Roger Stafford

Date: 14 Dec, 2011 05:17:07

Message: 3 of 5

"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

Subject: Finite difference method (implicit) for 1-D heat equation

From: Roger Stafford

Date: 14 Dec, 2011 16:36:08

Message: 4 of 5

"Roger Stafford" wrote in message <jc9bgj$8b5$1@newscl01ah.mathworks.com>...
> 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.) ......
- - - - - - - - - -
  After a good night's sleep I see that I misled you on the finite differences I described. I should have been dividing by h^2 and k rather than multiplying. I'm sorry about that. It was a silly mistake.

  Also apparently in your computation you are using U(i,~) rather than U(i-1,~) in computing the second difference, which would give you these two quantities:

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

and

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

to be equated. That would also make far more sense than the version I gave which would totally ignore changes in time occurring at the two x-boundaries.

  However, it seems to me that greater accuracy ought to be obtained by the average of these two second differences since (U(i,j)-U(i-1,j))/k most accurately estimates the partial derivative, dU/dt, at the midpoint between i and i-1 rather than at i. In other words in my opinion these two quantities should be equal:

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

and

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

  I don't believe that is what your procedure carries out. You might consider trying this to see if the accuracy improves. It would of course require a revision of your M and P matrices' definitions.

Roger Stafford

Subject: Finite difference method (implicit) for 1-D heat equation

From: Demoslav

Date: 14 Dec, 2011 17:01:08

Message: 5 of 5

"Roger Stafford" wrote in message <jcaj9o$b7m$1@newscl01ah.mathworks.com>...
> "Roger Stafford" wrote in message <jc9bgj$8b5$1@newscl01ah.mathworks.com>...
> > 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.) ......
> - - - - - - - - - -
> After a good night's sleep I see that I misled you on the finite differences I described. I should have been dividing by h^2 and k rather than multiplying. I'm sorry about that. It was a silly mistake.
>
> Also apparently in your computation you are using U(i,~) rather than U(i-1,~) in computing the second difference, which would give you these two quantities:
>
> (U(i,j-1)-2*U(i,j)+U(i,j+1))/h^2
>
> and
>
> a2*(U(i,j)-U(i-1,j))/k
>
> to be equated. That would also make far more sense than the version I gave which would totally ignore changes in time occurring at the two x-boundaries.
>
> However, it seems to me that greater accuracy ought to be obtained by the average of these two second differences since (U(i,j)-U(i-1,j))/k most accurately estimates the partial derivative, dU/dt, at the midpoint between i and i-1 rather than at i. In other words in my opinion these two quantities should be equal:
>
> 1/2*(U(i-1,j-1)-2*U(i-1,j)+U(i-1,j+1)+U(i,j-1)-2*U(i,j)+U(i,j+1))/h^2
>
> and
>
> a2*(U(i,j)-U(i-1,j))/k
>
> I don't believe that is what your procedure carries out. You might consider trying this to see if the accuracy improves. It would of course require a revision of your M and P matrices' definitions.
>
> Roger Stafford

you´re right, but then it would be crank - nicholson and my task is the implicit method. I ve made a revision of that for cycle and i figured out that my code is OK. i computed the relative error ( or what ever it is :) by substracting the solutions, then i made the frobenius norm of that matrix divided by the frob norm of the num solution matrix (i hope u understand - my math english isnt very good). the valus were aroun 0,001-0,0001 so i supose it is ok. the point is that i think (i dont know wheter it is right) that when the solution has big values (eg 10^6) then the error can be around 100 (???).
...but my code should be OK:)

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us