Asked by Abed Alnaif
on 8 Dec 2012

Hi,

I've been trying to solve a non-linear, heat-equation-type system of PDE's using the 'pdepe' function, with only one dimension in space. However, for many sets of parameter values, the solver exhibits unstable behaviour (oscillations, etc). I think that my problem demands a more sophisticated solver, due to nonlinearities and discontinuities. Hence, I am now trying to use the PDE Toolbox, hoping that it would be able to handle my problem, since it has an adaptive mesh algorithm, etc. However, I was unable to figure out how to use the PDE Toolbox for one-dimensional problems, nor was I able to find any examples of this.

Is it possible to use the PDE Toolbox to solve one-dimensional problems? I'm find with using the command-line interface (I don't necessarily need to use the GUI). If this is possible, I'd greatly appreciate it if someone could provide me with some code that solves a simple heat equation PDE (one dimension in space) using the PDE Toolbox, just to get me on the right foot (since, for some reason, I get really lost when I try to read the documentation for the PDE Toolbox).

Thanks,

Abed

*No products are associated with this question.*

Answer by Bill Greene
on 8 Dec 2012

Accepted answer

Hi,

Solving a 1-D PDE with PDE Toolbox is fairly straightforward. You just define a rectangular region of the appropriate width and arbitrary height. On the two edges at y=constant, you want a zero-Neumann BC. On the edges at x=constant, the BCs should not vary in y.

I've appended a very simple example of time-dependent heat transfer in a bar below.

Based on your description of the problem you are trying to solve, however, I can't think of any reason why PDE Toolbox should give a better solution than pdepe.

Bill

function simple1DTest % 1D transient heat transfer in x-direction h =.1; w=1; % width equal one, height is arbitrary (set to .1) g = decsg([3 4 0 w w 0 0 0 h h]', 'R1', ('R1')'); [p, e, t]=initmesh(g, 'Hmax', .05); b=@boundFile; c=1; a=0; d=1; f='100*x'; % heat load varies along the bar u0 = 0; % initial temperature equals zero tlist = 0:.02:1; u=parabolic(u0, tlist, b,p,e,t,c,a,f,d); figure; pdeplot(p,e,t, 'xydata', u(:,end), 'contour', 'on'); axis equal; title(sprintf('Temperature Distribution at time=%g seconds', tlist(end))); figure; plot(tlist, u(2,:)); xlabel 'Time'; ylabel 'Temperature'; grid; title 'Temperature at tip as a function of time' end

function [ q, g, h, r ] = boundFile( p, e, u, time ) N = 1; ne = size(e,2); q = zeros(N^2, ne); g = zeros(N, ne); h = zeros(N^2, 2*ne); r = zeros(N, 2*ne); % zero Neumann BCs (insulated) on edges 1 and 3 at y=0 and y=h % and the right edge, edge 2 for i=1:ne switch(e(5,i)) case 4 h(1,i) = 1; h(1,i+ne) = 1; r(i) = 500; r(i+ne) = 500; % 500 on left edge end end end

Answer by tom_brot
on 20 Mar 2015

Edited by tom_brot
on 20 Mar 2015

Ok this question has already been answered, but I want to share my experience with pdepe.

I also had problems with oscilating, unstable behaviour of pdepe. I found that the Problem lies not directly in the pdepe, but in the underlying ode15s solver. You can specify some options (with odeset()) for the pdepe-solver which are then passed to the ode15s solver.

The problem is, that only very few of the possible odeset() options are passed. To get a much more stable solution from the underlying ode15s solver I reduced the 'MaxOrder' option from the Default value of {5} to 2.

It's not that hard to make the option available for your pdepe solver.

**How to make your pdepe solver more stable:**

1) Find the file 'pdepe.m' in somewhere in your matlab installation path and copy it to your working directory. And rename it to e.g. 'pdepe_opt.m'.

2) Find the file 'pdentrp.m' (which is located in the 'private' directory next to 'pdepe.m') and copy it to your working directory.

3) In your local copy of pdepe.m, namely pdepe_opt.m add the line

maxorder = odeget(options,'MaxOrder',[],'fast');

after the line

maxstep = odeget(options,'MaxStep',[],'fast');

4) Then add

'MaxOrder',maxorder,

as argument of the odeset command, which is some lines below.

5) Before the line in your code where you call pdepe, insert

options = odeset('MaxOrder',2)

(or any number from 1 to 5) and then call the optimized pdepe function with your specified options:

sol = pdepe_opt(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t,options)

I hope this works as well for you as it did for me.

Best regards,

Tom

Opportunities for recent engineering grads.

## 0 Comments