Partial Differential equation code help!

1 view (last 30 days)
Mo
Mo on 14 May 2014
Commented: Bill Greene on 14 May 2014
Trying to solve the following partial differential equation:
I have the code:
function pdex1
m = 0;
x = linspace(0,1,20);
t = linspace(0,2,5);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
u = sol(:,:,1);
% A surface plot is often a good way to study a solution.
surf(x,t,u)
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')
% A solution profile can also be illuminating.
figure
plot(x,u(end,:))
title('Solution at t = 2')
xlabel('Distance x')
ylabel('u(x,2)')
% --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
D=3e-09;
v=1e-10;
k=1e3;
c = 1;
f = -D*(DuDx);
s = v-(k*u^2);
% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = 6.60e-07;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;
The problem I am having is as follows:
Where the equation was re-arranged in the form given. But I do not know how to put the term:
v-k[NO]^2
Into the code. Any help would be appreciated.
  4 Comments
Star Strider
Star Strider on 14 May 2014
Those are simply warnings telling you that the solver is having problems.
As for the errors in the surf call, I repeat: ‘check the elements of sol and the size of u to be sure u is a matrix with an appropriate size for x, and t’. I suspect the reason for the surf problem is that the solver is stopping early and not solving the equation on the mesh. That may be due to the extremely small magnitudes of D and v. Multiply them (and k) by 1E+3 or 1E+6 and see if the problem persists.
You also copied pdex1bc from the example. It may not be appropriate to your equation.
Bill Greene
Bill Greene on 14 May 2014
I agree with Mr. Strider's comment about there being a problem with your boundary conditions. You need to replace pdex1bc with one appropriate to your problem.
To understand the reason why the solver can't continue, I suggest doing the following. If you replace your line 3 with something like this:
t = linspace(0, .02, 100);
Then call pdepe the same way as you do now and add these lines right after it:
nt = size(sol, 1);
figure; plot(t(1:nt), sol(:,end));
you will see from the plot that the solution at the right end of the domain is getting very large; that is the reason the ODE solver can't continue.
Bill

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!