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:
Partial Differential Heat Equation

Subject: Partial Differential Heat Equation

From: Scott

Date: 15 Aug, 2010 14:00:05

Message: 1 of 23

This is my question:

Use the finite difference appraoch with h=0.2 and k=0.02 to solve the following partial differential heat equation using both the explicit and implicit method:

Note: I'll substitute the symbol for partial derivatives with a 'd' for simplicity...

d^2u/dx^2 = du/dt

The variable u(x,t) is analogous to temperature.
The range of distance x is 0<x<2 and time span is 0<t<0.5
The boundary conditions are u(o,t)=2 and u(2,t)=2
The initial condition is u(x,0)=2(x-1)^3

What is a good starting point to do this on MATLab?

Subject: Partial Differential Heat Equation

From: Roland

Date: 15 Aug, 2010 15:48:03

Message: 2 of 23

"Scott " <boofheadfenn@hotmail.com> wrote in message <i48rt5$n4r$1@fred.mathworks.com>...
> This is my question:
>
> Use the finite difference appraoch with h=0.2 and k=0.02 to solve the following partial differential heat equation using both the explicit and implicit method:
>
> Note: I'll substitute the symbol for partial derivatives with a 'd' for simplicity...
>
> d^2u/dx^2 = du/dt
>
> The variable u(x,t) is analogous to temperature.
> The range of distance x is 0<x<2 and time span is 0<t<0.5
> The boundary conditions are u(o,t)=2 and u(2,t)=2
> The initial condition is u(x,0)=2(x-1)^3
>
> What is a good starting point to do this on MATLab?

Well, I think numerical methods are a topic on its own. We made a CFD code at University in Matlab, maybe i can give you a rough idea about what to do.

- discretice your partial differential heat equation on your own
- set up a mesh, in your case 1D, mesh means: a vector with x values
- set up a matrix for u, i would make it a 2D Matrix, with as many colums as you have cells in your mesh and number of rows are your time steps
- explicitely solved, loop through the cells any apply the descreticed heat equation for every cell, until cell values dont change much (converge)
- implicetely solved, you have to solve a set of linear equations

Subject: Partial Differential Heat Equation

From: Scott

Date: 19 Aug, 2010 07:38:08

Message: 3 of 23

"Roland " <burgmann@gmx.de> wrote in message <i4927j$71h$1@fred.mathworks.com>...
> "Scott " <boofheadfenn@hotmail.com> wrote in message <i48rt5$n4r$1@fred.mathworks.com>...
> > This is my question:
> >
> > Use the finite difference appraoch with h=0.2 and k=0.02 to solve the following partial differential heat equation using both the explicit and implicit method:
> >
> > Note: I'll substitute the symbol for partial derivatives with a 'd' for simplicity...
> >
> > d^2u/dx^2 = du/dt
> >
> > The variable u(x,t) is analogous to temperature.
> > The range of distance x is 0<x<2 and time span is 0<t<0.5
> > The boundary conditions are u(o,t)=2 and u(2,t)=2
> > The initial condition is u(x,0)=2(x-1)^3
> >
> > What is a good starting point to do this on MATLab?
>
> Well, I think numerical methods are a topic on its own. We made a CFD code at University in Matlab, maybe i can give you a rough idea about what to do.
>
> - discretice your partial differential heat equation on your own
> - set up a mesh, in your case 1D, mesh means: a vector with x values
> - set up a matrix for u, i would make it a 2D Matrix, with as many colums as you have cells in your mesh and number of rows are your time steps
> - explicitely solved, loop through the cells any apply the descreticed heat equation for every cell, until cell values dont change much (converge)
> - implicetely solved, you have to solve a set of linear equations

Thanks Roland,

I sort of get what you're saying. I'm new to MATLab so I'm sorry about this. I tried what I think is the Explicit Method and I'll post my coding below. I have compared the surface plot to a surface plot I created using Excel and it is not right. Can you see where I might be going wrong? I'm pretty sure it has something to do with my boundary conditions because on the surface plot my boundary condition for length (x)=0 is 2 so all my u values should be 2 but they aren't, and likewise for when length (x)=2.

function pdex1

m = 0;
x = linspace(0,2,11);
t = linspace(0,0.5,26);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
u = sol(:,:,1);

figure;
surf(x,t,u);
title('Numerical solution');
xlabel('Distance x');
ylabel('Time t');

function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = 1;
f = DuDx;
s = 0;

function u0 = pdex1ic(x)
u0 = 2*((x-1)^3);

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur;
qr = 2; (Pretty sure there is something wrong with these values??)

Subject: Partial Differential Heat Equation

From: Torsten Hennig

Date: 19 Aug, 2010 08:04:06

Message: 4 of 23

> "Roland " <burgmann@gmx.de> wrote in message
> <i4927j$71h$1@fred.mathworks.com>...
> > "Scott " <boofheadfenn@hotmail.com> wrote in
> message <i48rt5$n4r$1@fred.mathworks.com>...
> > > This is my question:
> > >
> > > Use the finite difference appraoch with h=0.2 and
> k=0.02 to solve the following partial differential
> heat equation using both the explicit and implicit
> method:
> > >
> > > Note: I'll substitute the symbol for partial
> derivatives with a 'd' for simplicity...
> > >
> > > d^2u/dx^2 = du/dt
> > >
> > > The variable u(x,t) is analogous to temperature.
> > > The range of distance x is 0<x<2 and time span is
> 0<t<0.5
> > > The boundary conditions are u(o,t)=2 and u(2,t)=2
> > > The initial condition is u(x,0)=2(x-1)^3
> > >
> > > What is a good starting point to do this on
> MATLab?
> >
> > Well, I think numerical methods are a topic on its
> own. We made a CFD code at University in Matlab,
> maybe i can give you a rough idea about what to do.
> >
> > - discretice your partial differential heat
> equation on your own
> > - set up a mesh, in your case 1D, mesh means: a
> vector with x values
> > - set up a matrix for u, i would make it a 2D
> Matrix, with as many colums as you have cells in your
> mesh and number of rows are your time steps
> > - explicitely solved, loop through the cells any
> apply the descreticed heat equation for every cell,
> until cell values dont change much (converge)
> > - implicetely solved, you have to solve a set of
> linear equations
>
> Thanks Roland,
>
> I sort of get what you're saying. I'm new to MATLab
> so I'm sorry about this. I tried what I think is the
> Explicit Method and I'll post my coding below. I have
> compared the surface plot to a surface plot I created
> using Excel and it is not right. Can you see where I
> might be going wrong? I'm pretty sure it has
> something to do with my boundary conditions because
> on the surface plot my boundary condition for length
> (x)=0 is 2 so all my u values should be 2 but they
> aren't, and likewise for when length (x)=2.
>
> function pdex1
>
> m = 0;
> x = linspace(0,2,11);
> t = linspace(0,0.5,26);
> sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
> u = sol(:,:,1);
>
> figure;
> surf(x,t,u);
> title('Numerical solution');
> xlabel('Distance x');
> ylabel('Time t');
>
> function [c,f,s] = pdex1pde(x,t,u,DuDx)
> c = 1;
> f = DuDx;
> s = 0;
>
> function u0 = pdex1ic(x)
> u0 = 2*((x-1)^3);
>
> function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
> pl = ul;
> ql = 0;
> pr = ur;
> qr = 2; (Pretty sure there is something wrong
> with these values??)

At t=0, your initial condition gives u=-2 at x=0.
Is it that what you want ?
Then the initial condition is inconsistent with
the boundary condition (which will result in
difficulties for the solver).

The boundary conditions routine for fixed values
of u should look like
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul-ulb;
ql = 0;
pr = ur-urb;
qr = 2
where ulb and urb are the values you want to prescribe
for u at the left resp. right boundary point.

Best wishes
Torsten.

Subject: Partial Differential Heat Equation

From: Torsten Hennig

Date: 19 Aug, 2010 08:17:40

Message: 5 of 23

> > "Roland " <burgmann@gmx.de> wrote in message
> > <i4927j$71h$1@fred.mathworks.com>...
> > > "Scott " <boofheadfenn@hotmail.com> wrote in
> > message <i48rt5$n4r$1@fred.mathworks.com>...
> > > > This is my question:
> > > >
> > > > Use the finite difference appraoch with h=0.2
> and
> > k=0.02 to solve the following partial differential
> > heat equation using both the explicit and implicit
> > method:
> > > >
> > > > Note: I'll substitute the symbol for partial
> > derivatives with a 'd' for simplicity...
> > > >
> > > > d^2u/dx^2 = du/dt
> > > >
> > > > The variable u(x,t) is analogous to
> temperature.
> > > > The range of distance x is 0<x<2 and time span
> is
> > 0<t<0.5
> > > > The boundary conditions are u(o,t)=2 and
> u(2,t)=2
> > > > The initial condition is u(x,0)=2(x-1)^3
> > > >
> > > > What is a good starting point to do this on
> > MATLab?
> > >
> > > Well, I think numerical methods are a topic on
> its
> > own. We made a CFD code at University in Matlab,
> > maybe i can give you a rough idea about what to
> do.
> > >
> > > - discretice your partial differential heat
> > equation on your own
> > > - set up a mesh, in your case 1D, mesh means: a
> > vector with x values
> > > - set up a matrix for u, i would make it a 2D
> > Matrix, with as many colums as you have cells in
> your
> > mesh and number of rows are your time steps
> > > - explicitely solved, loop through the cells any
> > apply the descreticed heat equation for every
> cell,
> > until cell values dont change much (converge)
> > > - implicetely solved, you have to solve a set of
> > linear equations
> >
> > Thanks Roland,
> >
> > I sort of get what you're saying. I'm new to
> MATLab
> > so I'm sorry about this. I tried what I think is
> the
> > Explicit Method and I'll post my coding below. I
> have
> > compared the surface plot to a surface plot I
> created
> > using Excel and it is not right. Can you see where
> I
> > might be going wrong? I'm pretty sure it has
> > something to do with my boundary conditions
> because
> > on the surface plot my boundary condition for
> length
> > (x)=0 is 2 so all my u values should be 2 but they
> > aren't, and likewise for when length (x)=2.
> >
> > function pdex1
> >
> > m = 0;
> > x = linspace(0,2,11);
> > t = linspace(0,0.5,26);
> > sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
> > u = sol(:,:,1);
> >
> > figure;
> > surf(x,t,u);
> > title('Numerical solution');
> > xlabel('Distance x');
> > ylabel('Time t');
> >
> > function [c,f,s] = pdex1pde(x,t,u,DuDx)
> > c = 1;
> > f = DuDx;
> > s = 0;
> >
> > function u0 = pdex1ic(x)
> > u0 = 2*((x-1)^3);
> >
> > function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
> > pl = ul;
> > ql = 0;
> > pr = ur;
> > qr = 2; (Pretty sure there is something
> wrong
> > with these values??)
>
> At t=0, your initial condition gives u=-2 at x=0.
> Is it that what you want ?
> Then the initial condition is inconsistent with
> the boundary condition (which will result in
> difficulties for the solver).
>
> The boundary conditions routine for fixed values
> of u should look like
> function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
> pl = ul-ulb;
> ql = 0;
> pr = ur-urb;
> qr = 2
      
Sorry, should be
qr=0;
 
> where ulb and urb are the values you want to
> prescribe
> for u at the left resp. right boundary point.
>
> Best wishes
> Torsten.

Subject: Partial Differential Heat Equation

From: Scott

Date: 19 Aug, 2010 08:40:20

Message: 6 of 23

Hi Torsten,

No my Boundary Conditions are:

u(0,t)=2
u(2,t)=2

and my Initial Condition is:

u(x,0)=2(x-1)^3

Does this change your answer much?

Regards Scott

Subject: Partial Differential Heat Equation

From: Scott

Date: 19 Aug, 2010 08:45:09

Message: 7 of 23

So I have,

∂^2u/∂x^2 = ∂u/∂t

h=0.2, k=0.02
0<x<2 and 0<t<0.5
with the previously mentioned conditions.

I have to solve using explicit method and implicit method

Regards,
Scott

Subject: Partial Differential Heat Equation

From: Scott

Date: 19 Aug, 2010 09:39:05

Message: 8 of 23

Also, any chance you know what else might be wrong with my coding? When I run it the following errors come up:

??? Undefined function or variable 'ulb'.

Error in ==> pdex1>pdex1bc at 86
pl = ul-ulb;

Error in ==> pdepe at 266
[pL,qL,pR,qR] = feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),t(1),varargin{:});

Error in ==> pdex1 at 45
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

Error in ==> run at 57
          evalin('caller', [s ';']);

Subject: Partial Differential Heat Equation

From: Torsten Hennig

Date: 19 Aug, 2010 10:02:22

Message: 9 of 23

> Also, any chance you know what else might be wrong
> with my coding? When I run it the following errors
> come up:
>
> ??? Undefined function or variable 'ulb'.
>
> Error in ==> pdex1>pdex1bc at 86
> pl = ul-ulb;
>
> Error in ==> pdepe at 266
> [pL,qL,pR,qR] =
> feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),t(1),vara
> rgin{:});
>
> Error in ==> pdex1 at 45
> sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
>
> Error in ==> run at 57
> evalin('caller', [s ';']);

If your boundary conditions are u=2 at both ends,
just substitute 2 for ulr and ulb.

If your initial condition is really u(x)=2*(x-1)^3,
prescribe this condition for all x except x=0
and prescribe u=2 at x=0 in the pdexic-routine.

Best wishes
Torsten.

Subject: Partial Differential Heat Equation

From: Torsten Hennig

Date: 19 Aug, 2010 10:07:31

Message: 10 of 23

> > Also, any chance you know what else might be wrong
> > with my coding? When I run it the following errors
> > come up:
> >
> > ??? Undefined function or variable 'ulb'.
> >
> > Error in ==> pdex1>pdex1bc at 86
> > pl = ul-ulb;
> >
> > Error in ==> pdepe at 266
> > [pL,qL,pR,qR] =
> >
> feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),t(1),vara
>
> > rgin{:});
> >
> > Error in ==> pdex1 at 45
> > sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
> >
> > Error in ==> run at 57
> > evalin('caller', [s ';']);
>
> If your boundary conditions are u=2 at both ends,
> just substitute 2 for ulr and ulb.
>

 If your boundary conditions are u=2 at both ends,
 just substitute 2 for urb and ulb.

> If your initial condition is really u(x)=2*(x-1)^3,
> prescribe this condition for all x except x=0
> and prescribe u=2 at x=0 in the pdexic-routine.
>
> Best wishes
> Torsten.

Best wishes
Torsten.

Subject: Partial Differential Heat Equation

From: Scott

Date: 19 Aug, 2010 10:22:06

Message: 11 of 23

Wow, thanks so much for that... I can't believe it was that easy!

Do you have any idea how this could be solved using the implicit method?

Subject: Partial Differential Heat Equation

From: Scott

Date: 19 Aug, 2010 10:35:04

Message: 12 of 23

Wow, thanks so much for that... I can't believe it was that easy!

Do you have any idea how this could be solved using the implicit method?

Subject: Partial Differential Heat Equation

From: Torsten Hennig

Date: 19 Aug, 2010 10:43:00

Message: 13 of 23

> Wow, thanks so much for that... I can't believe it
> was that easy!
>
> Do you have any idea how this could be solved using
> the implicit method?


I think what you did is not what you are expected to
do.

Take a look at
http://en.wikipedia.org/wiki/Finite_difference_method
subsection
Example: The heat equation
Explicit method
Implicit method

to see how both methods work.

pdepe in MATLAB is a good tool to check your results
obtained by the two methods described above, but
it is neither the explicit nor the implicit method ;
it's the method of lines.

Best wishes
Torsten.

Subject: Partial Differential Heat Equation

From: Scott

Date: 19 Aug, 2010 11:07:04

Message: 14 of 23

Hmm as much as I'd hate to admit it I think you're right.

Is there any advice you can give me how to solve the equation using the explicit or implicit methods? Or any examples I could follow?

Many thanks,
Scott

Subject: Partial Differential Heat Equation

From: Scott

Date: 19 Aug, 2010 12:49:06

Message: 15 of 23

Ok, So I've had a go at the explicit version and the results I'm getting are exactly the same as when I did my Excel model up which is good, however it is taking me forever to do and I'm not sure how I'd ever be able to graph it because of all my different u values. I'll show you my code, can you tell me if there's any way I can shorten this procedure so that I don't have to repeat and repeat and repeat. I think I need to do a for loop of some sort which I'm very rusty on. Anyway here's my code:

h=0.2;
k=0.02;
% range of distance x is 0<x<2, divided into proportions of h:
x=(0: h: 2);
% range of time span is 0<t<0.5 divided into proportions of k:
t=(0: k: 0.5);

% boundary condition u(0,t)=2:
u0_0=2;
u0_1=2;
u0_2=2;
u0_3=2;
u0_4=2;
u0_5=2;
u0_6=2;
u0_7=2;
u0_8=2;
u0_9=2;
u0_10=2;
u0_11=2;
u0_12=2;
u0_13=2;
u0_14=2;
u0_15=2;
u0_16=2;
u0_17=2;
u0_18=2;
u0_19=2;
u0_20=2;
u0_21=2;
u0_22=2;
u0_23=2;
u0_24=2;
u0_25=2;

% boundary condition u(2,t)=2:
u10_0=2;
u10_1=2;
u10_2=2;
u10_3=2;
u10_4=2;
u10_5=2;
u10_6=2;
u10_7=2;
u10_8=2;
u10_9=2;
u10_10=2;
u10_11=2;
u10_12=2;
u10_13=2;
u10_14=2;
u10_15=2;
u10_16=2;
u10_17=2;
u10_18=2;
u10_19=2;
u10_20=2;
u10_21=2;
u10_22=2;
u10_23=2;
u10_24=2;
u10_25=2;

% initial condition u(x,0) = 2*(x-1)^3 for t=0 and 0<=x<=2:

u1_0=2*(x(2)-1)^3;
u2_0=2*(x(3)-1)^3;
u3_0=2*(x(4)-1)^3;
u4_0=2*(x(5)-1)^3;
u5_0=2*(x(6)-1)^3;
u6_0=2*(x(7)-1)^3;
u7_0=2*(x(8)-1)^3;
u8_0=2*(x(9)-1)^3;
u9_0=2*(x(10)-1)^3;

% all intermediate u values following the formula:
% u(j,n+1)=(1-2*lambda)u(j,n)+lambda*u(j-1,n)+lambda*u(j+1,n)

lambda=k/(h^2);

% for t=(2):

u1_1=(1-2*lambda)*u1_0+lambda*u0_0+lambda*u2_0;
u2_1=(1-2*lambda)*u2_0+lambda*u1_0+lambda*u3_0;
u3_1=(1-2*lambda)*u3_0+lambda*u2_0+lambda*u4_0;
u4_1=(1-2*lambda)*u4_0+lambda*u3_0+lambda*u5_0;
u5_1=(1-2*lambda)*u5_0+lambda*u4_0+lambda*u6_0;
u6_1=(1-2*lambda)*u6_0+lambda*u5_0+lambda*u7_0;
u7_1=(1-2*lambda)*u7_0+lambda*u6_0+lambda*u8_0;
u8_1=(1-2*lambda)*u8_0+lambda*u7_0+lambda*u9_0;
u9_1=(1-2*lambda)*u9_0+lambda*u8_0+lambda*u10_0;

% for t=(3):

u1_2=(1-2*lambda)*u1_1+lambda*u0_1+lambda*u2_1;
u2_2=(1-2*lambda)*u2_1+lambda*u1_1+lambda*u3_1;
u3_2=(1-2*lambda)*u3_1+lambda*u2_1+lambda*u4_1;
u4_2=(1-2*lambda)*u4_1+lambda*u3_1+lambda*u5_1;
u5_2=(1-2*lambda)*u5_1+lambda*u4_1+lambda*u6_1;
u6_2=(1-2*lambda)*u6_1+lambda*u5_1+lambda*u7_1;
u7_2=(1-2*lambda)*u7_1+lambda*u6_1+lambda*u8_1;
u8_2=(1-2*lambda)*u8_1+lambda*u7_1+lambda*u9_1;
u9_2=(1-2*lambda)*u9_1+lambda*u8_1+lambda*u10_1;

etc etc etc

Subject: Partial Differential Heat Equation

From: Torsten Hennig

Date: 19 Aug, 2010 13:37:55

Message: 16 of 23

> Ok, So I've had a go at the explicit version and the
> results I'm getting are exactly the same as when I
> did my Excel model up which is good, however it is
> taking me forever to do and I'm not sure how I'd ever
> be able to graph it because of all my different u
> values. I'll show you my code, can you tell me if
> there's any way I can shorten this procedure so that
> I don't have to repeat and repeat and repeat. I think
> I need to do a for loop of some sort which I'm very
> rusty on. Anyway here's my code:
>
> h=0.2;
> k=0.02;
> % range of distance x is 0<x<2, divided into
> proportions of h:
> x=(0: h: 2);
> % range of time span is 0<t<0.5 divided into
> proportions of k:
> t=(0: k: 0.5);
>
> % boundary condition u(0,t)=2:
> u0_0=2;
> u0_1=2;
> u0_2=2;
> u0_3=2;
> u0_4=2;
> u0_5=2;
> u0_6=2;
> u0_7=2;
> u0_8=2;
> u0_9=2;
> u0_10=2;
> u0_11=2;
> u0_12=2;
> u0_13=2;
> u0_14=2;
> u0_15=2;
> u0_16=2;
> u0_17=2;
> u0_18=2;
> u0_19=2;
> u0_20=2;
> u0_21=2;
> u0_22=2;
> u0_23=2;
> u0_24=2;
> u0_25=2;
>
> % boundary condition u(2,t)=2:
> u10_0=2;
> u10_1=2;
> u10_2=2;
> u10_3=2;
> u10_4=2;
> u10_5=2;
> u10_6=2;
> u10_7=2;
> u10_8=2;
> u10_9=2;
> u10_10=2;
> u10_11=2;
> u10_12=2;
> u10_13=2;
> u10_14=2;
> u10_15=2;
> u10_16=2;
> u10_17=2;
> u10_18=2;
> u10_19=2;
> u10_20=2;
> u10_21=2;
> u10_22=2;
> u10_23=2;
> u10_24=2;
> u10_25=2;
>
> % initial condition u(x,0) = 2*(x-1)^3 for t=0 and
> 0<=x<=2:
>
> u1_0=2*(x(2)-1)^3;
> u2_0=2*(x(3)-1)^3;
> u3_0=2*(x(4)-1)^3;
> u4_0=2*(x(5)-1)^3;
> u5_0=2*(x(6)-1)^3;
> u6_0=2*(x(7)-1)^3;
> u7_0=2*(x(8)-1)^3;
> u8_0=2*(x(9)-1)^3;
> u9_0=2*(x(10)-1)^3;
>
> % all intermediate u values following the formula:
> %
> u(j,n+1)=(1-2*lambda)u(j,n)+lambda*u(j-1,n)+lambda*u(j
> +1,n)
>
> lambda=k/(h^2);
>
> % for t=(2):
>
> u1_1=(1-2*lambda)*u1_0+lambda*u0_0+lambda*u2_0;
> u2_1=(1-2*lambda)*u2_0+lambda*u1_0+lambda*u3_0;
> u3_1=(1-2*lambda)*u3_0+lambda*u2_0+lambda*u4_0;
> u4_1=(1-2*lambda)*u4_0+lambda*u3_0+lambda*u5_0;
> u5_1=(1-2*lambda)*u5_0+lambda*u4_0+lambda*u6_0;
> u6_1=(1-2*lambda)*u6_0+lambda*u5_0+lambda*u7_0;
> u7_1=(1-2*lambda)*u7_0+lambda*u6_0+lambda*u8_0;
> u8_1=(1-2*lambda)*u8_0+lambda*u7_0+lambda*u9_0;
> u9_1=(1-2*lambda)*u9_0+lambda*u8_0+lambda*u10_0;
>
> % for t=(3):
>
> u1_2=(1-2*lambda)*u1_1+lambda*u0_1+lambda*u2_1;
> u2_2=(1-2*lambda)*u2_1+lambda*u1_1+lambda*u3_1;
> u3_2=(1-2*lambda)*u3_1+lambda*u2_1+lambda*u4_1;
> u4_2=(1-2*lambda)*u4_1+lambda*u3_1+lambda*u5_1;
> u5_2=(1-2*lambda)*u5_1+lambda*u4_1+lambda*u6_1;
> u6_2=(1-2*lambda)*u6_1+lambda*u5_1+lambda*u7_1;
> u7_2=(1-2*lambda)*u7_1+lambda*u6_1+lambda*u8_1;
> u8_2=(1-2*lambda)*u8_1+lambda*u7_1+lambda*u9_1;
> u9_2=(1-2*lambda)*u9_1+lambda*u8_1+lambda*u10_1;
>
> etc etc etc


I won't write the code for you, but you should try
to become familiar with arrays (use a two-dimensional array for u (first index: counter for timestep,
second index: counter for position)) and for-loops
(use a double loop (outer loop: time, inner loop: position)).
 
Best wishes
Torsten.

Subject: Partial Differential Heat Equation

From: Scott

Date: 19 Aug, 2010 13:53:26

Message: 17 of 23

Ok, from the wikipedia example I've got this code now for Explicit method:

function heat(f,c1,c2,L,T,h,k,alpha)
n=L/h;
m=T/k;
lambda=alpha^2*k/(h^2)
x=0:h:L;
for i=1:n+1
    u(i)=feval(f,(i-1)*h);
end
for j=1:m
    t=j*k;
    for i=1:n+1
        if (i==1)
            y(i)=c1;
        elseif (i==n+1)
            y(i)=c2;
        else
            y(i)=(1-2*lambda)*u(i)+lambda*(u(i+1)+u(i-1));
        end;
    end;
    u=y;
end


and in a separate m file i have this:

f=2*(x-1)^3;
c1=2;
c2=2;
L=2;
T=0.5;
h=0.2;
k=0.02;
alpha=1;

I can't seem to run either without there being an error. Where am I going wrong???
This is killing me!!

Subject: Partial Differential Heat Equation

From: Roland

Date: 20 Aug, 2010 13:50:09

Message: 18 of 23

this at least does to anything:

clc
clear all

f=@(x)2.*(x-1).^3;
c1=2;
c2=2;
L=2;
T=0.5;
h=0.2;
k=0.02;
alpha=1;
n=L/h;
m=T/k;

lambda=alpha^2*k/(h^2);
x=0:h:L;
for i=1:n+1
    u(i)=feval(f,(i-1)*h);
end
for j=1:m
    t=j*k;
    for i=1:n+1
        if (i==1)
            y(i)=c1;
        elseif (i==n+1)
            y(i)=c2;
        else
            y(i)=(1-2*lambda)*u(i)+lambda*(u(i+1)+u(i-1));
        end;
    end;
    u=y;
end

Subject: Partial Differential Heat Equation

From: Jan-Pou Nee

Date: 11 Feb, 2011 05:27:03

Message: 19 of 23

"SF88" wrote in message <i4iqlk$esc$1@fred.mathworks.com>...
> Hi Torsten,
>
> No my Boundary Conditions are:
>
> u(0,t)=2
> u(2,t)=2
>
> and my Initial Condition is:
>
> u(x,0)=2(x-1)^3
>
> Does this change your answer much?
>
> Regards Scott
What if my initial value is a step function
u(x,0)=5 if 0<=x<0.5 and
u(x,0)=1 if 0.5<=x<=1
I use iu =@(x) (x <= 0.5).*5 + (0.5 <= x <= 1).*0;
the result doesn't looks like a step function
Then I create a vector iu but I don't know how to
plug this vector in using pdex1ic

Subject: Partial Differential Heat Equation

From: Torsten Hennig

Date: 11 Feb, 2011 07:23:23

Message: 20 of 23

> What if my initial value is a step function
> u(x,0)=5 if 0<=x<0.5 and
> u(x,0)=1 if 0.5<=x<=1
> I use iu =@(x) (x <= 0.5).*5 + (0.5 <= x <= 1).*0;
> the result doesn't looks like a step function
> Then I create a vector iu but I don't know how to
> plug this vector in using pdex1ic

Pass a function handle @iu to pdepe and define the
function iu as

function u0 = iu(x)
u0 = 0.0;
if (x<=0.5)
  u0 = 0.5;
end

Best wishes
Torsten.

Subject: Partial Differential Heat Equation

From: Torsten Hennig

Date: 11 Feb, 2011 09:35:45

Message: 21 of 23

> > What if my initial value is a step function
> > u(x,0)=5 if 0<=x<0.5 and
> > u(x,0)=1 if 0.5<=x<=1
> > I use iu =@(x) (x <= 0.5).*5 + (0.5 <= x <=
> 1).*0;
> > the result doesn't looks like a step function
> > Then I create a vector iu but I don't know how to
> > plug this vector in using pdex1ic
>
> Pass a function handle @iu to pdepe and define the
> function iu as
>
> function u0 = iu(x)
> u0 = 0.0;
> if (x<=0.5)
> u0 = 0.5;

Should read : u0 = 5;

> end
>
> Best wishes
> Torsten.

Subject: Partial Differential Heat Equation

From: Jan-Pou Nee

Date: 12 Feb, 2011 02:40:05

Message: 22 of 23

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <916009494.76673.1297416976196.JavaMail.root@gallium.mathforum.org>...
> > > What if my initial value is a step function
> > > u(x,0)=5 if 0<=x<0.5 and
> > > u(x,0)=1 if 0.5<=x<=1
> > > I use iu =@(x) (x <= 0.5).*5 + (0.5 <= x <=
> > 1).*0;
> > > the result doesn't looks like a step function
> > > Then I create a vector iu but I don't know how to
> > > plug this vector in using pdex1ic
> >
> > Pass a function handle @iu to pdepe and define the
> > function iu as
> >
> > function u0 = iu(x)
> > u0 = 0.0;
> > if (x<=0.5)
> > u0 = 0.5;
>
> Should read : u0 = 5;
>
> > end
> >
> > Best wishes
> > Torsten.
Thanks, it seems that the trial version has bugs on "pdex_"
I get solution with sin(x) as initial value, I am not sure though.

Subject: Partial Differential Heat Equation

From: Jan-Pou Nee

Date: 12 Feb, 2011 07:01:04

Message: 23 of 23

"Jan-Pou Nee" <jpntw@pchome.com.tw> wrote in message <ij4ru5$4cf$1@fred.mathworks.com>...
> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <916009494.76673.1297416976196.JavaMail.root@gallium.mathforum.org>...
> > > > What if my initial value is a step function
> > > > u(x,0)=5 if 0<=x<0.5 and
> > > > u(x,0)=1 if 0.5<=x<=1
> > > > I use iu =@(x) (x <= 0.5).*5 + (0.5 <= x <=
> > > 1).*0;
> > > > the result doesn't looks like a step function
> > > > Then I create a vector iu but I don't know how to
> > > > plug this vector in using pdex1ic
> > >
> > > Pass a function handle @iu to pdepe and define the
> > > function iu as
> > >
> > > function u0 = iu(x)
> > > u0 = 0.0;
> > > if (x<=0.5)
> > > u0 = 0.5;
> >
> > Should read : u0 = 5;
> >
> > > end
> > >
> > > Best wishes
> > > Torsten.
> Thanks, it seems that the trial version has bugs on "pdex_"
> I get solution with sin(x) as initial value, I am not sure though.
My mistake, set u0=pdex1ic(x) will do. Thanks again.

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