PDE propagating from point source

I want to solve numerically a nonlinear diffusion equation from an instantaneous point source. Thus, I have initial conditions, but not boundary conditions. How should I go about writing a code to solve circular propagation from a point?

Thanks!!

[Torsten]

What is the equation you try to solve (because you are talking about a nonlinear diffusion equation) ?

Best wishes

Torsten

[María Jesús]

$\frac{\partial C}{\partial t}=r^{1-s}\frac{\partial}{\partial r}[r^{s-1}D\frac{\partial C}{\partial r})]$ where $s$ is constant and $D=D_0(\frac{\partial C}{\partial C_0})^n$ and $n>0$

Answer 1

[Torsten]

I assume you want the point source appear at r=0.

Choose the interval of integration as [0:R] where R is big enough to ensure that C=C(t=0) throughout the period of integration.

As initial condition, choose an approximation to the delta function.

As boundary conditions, choose dC/dr = 0 at both ends.

Best wishes

Torsten.

[Torsten]

... and try MATLAB's pdepe to solve.

Best wishes

Torsten.

[María Jesús]

thanks! just one question: how do I specify the delta function? Also, do you mean Kronecker or dirac delta?

[Torsten]

Take a look at the paragraph "Approximations to the identity" under

https://en.wikipedia.org/wiki/Dirac_delta_function

Best wishes

Torsten.

[María Jesús]

I've written the code, but it gives errors... can anyone help with this? The code is

clear all

m = 1;

r = linspace(0, 1, 20);

t = linspace(0, 2, 10);

u = pdepe(m, @pat2, @ic1, @bc1, r, t);

surf(r, t, u);

title('Surface plot of solution');

xlabel('Distance r');

ylabel('Time t');

figure

plot(r, u(end,:))

title('Solution at t = 2')

xlabel('Distance r')

ylabel('u(r,2)')

function [c, b, s, D] = pat2(r, t, u, DuDr)

n = 1;

D_0 = 1;

D = D_0/(KronD(r, 0))^n;

c = 1;

b = D*u^n*DuDr;

s = 0;

function [d] = KronD(j,k)

if nargin < 2

    error('Too few inputs.  See help KronD')

elseif nargin > 2

    error('Too many inputs.  See help KronD')

end

if j == k

    d = Inf;

else

    d = 0;

end

function [pl, ql, pr, qr] = bc1(rl, ul, rr, ur, t)

pl = ul;

ql = 0;

pr = ur;

qr = 0;

function value = ic1(r)

value = KronD(r, 0);

And the error is

Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t. > In ode15s (line 668) In pdepe (line 289) In pattle2_0 (line 7) Warning: Time integration has failed. Solution is available at requested time points up to t=0.000000e+00. > In pdepe (line 303) In pattle2_0 (line 7) Error using surf (line 57) Z must be a matrix, not a scalar or vector.

Error in pattle2_0 (line 8) surf(r, t, u);

Thanks again!!!

[Torsten]

1. You will have to work with a numerical approximation to the delta function. I gave you a suitable link.

2. Your boundary conditions are incorrect. You will have to set

pl=0, ql=1, pr=0, qr=1

3. I don't understand your definition of D. The setting

D = D_0/(KronD(r, 0))^n;

doesn't make sense.

Best wishes

Torsten.

[Nicholas Mikolajewicz]

Torsten, regarding the earlier answer you provided, whats the reasoning behind using the dirac delta approximation for the point source rather than just setting the initial condition to the source concentration/density as u0(x==0) = initial condition?