How to solve this code error? (Diffusion_equation with pedpe)

% define the parameters
    D1=1;
    D2=2;
    a = 10;

   %define spatial mesh and time span
   x = linspace(0, a, 100);
   t = linspace(0, 10, 100);

   %use pdepe solver
   sol = pdepe(0, @diffusion_equation, @diffusion_ini, @diffusion_boun, x, t);

Unrecognized function or variable 'D1'.

Error in solution>diffusion_equation (line 28)
D(x <= 2/a) = D1; % Set D1 where the condition is true

Error in pdepe (line 242)
[c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});

   % Extract the solution
   u = sol(:,:,1);

   %plot the results
   figure;
   surf(x,t, u);
   xlabel('spatial coordinate (x)');
   ylabel('Time(t)');
   zlabel('Transport property profile (u)');
   title('1D diffusion model of the T cell density');


function [c, f, s] = diffusion_equation(x, t, u, DuDx)

    D = zeros(size(x));  % Initialize D as a vector of zeros
    D(x <= 2/a) = D1;     % Set D1 where the condition is true
    D(x > 2/a) = D2;      % Set D2 where the condition is false

    c = 1;
    f = D .* DuDx;
    s = 0;
end
%define the initial condition
function u0 = diffusion_ini(x)
    u0 = 10^3 * ones(size(x));
end
%define the boundary condition
function [pl,ql,pr,qr] = diffusion_boun(xl,ul,xr,ur,t)
    pl = ul;
    ql = 0;
    pr = ur;
    qr = 0;
end

this code is the diffusion equation but there is a condition of x

but there is code error. so i cannot run this code yet.

please help me

[Dyuman Joshi]

You're setting the elements to D1 and D2 here -

    D(x <= 2/a) = D1;     % Set D1 where the condition is true
    D(x > 2/a) = D2;      % Set D2 where the condition is false

But what are the values of D1 and D2?

[채현]

D1 and D2 are the diffusion coefficients in the diffusion equation, and what I want to see through this modeling is how the U(density) that I want to find is distributed as this diffusion coefficient, D, varies. So if you look at the x-axis, if the length is a, I arbitrarily set the coefficients to be D1 up to half of the length, 2/a, and D2 for the rest of the length, so D1<D2. At the top of the code, we've arbitrarily substituted the values 1 and 2.

The expected result is that D2 has a larger coefficient, so I'm expecting the graph to be denser towards this D2.

[Dyuman Joshi]

"At the top of the code, we've arbitrarily substituted the values 1 and 2."

Okay, yes. I seem to have overlooked those definitions.

Answer 1

[Dyuman Joshi]

Modify the function to provide D1, D2 and a as inputs and update the pdepe() call accordingly as well -

% define the parameters
    D1=1;
    D2=2;
    a = 10;

   %define spatial mesh and time span
   x = linspace(0, a, 100);
   t = linspace(0, 10, 100);

   %use pdepe solver
   %% Update the call to pdepde
   sol = pdepe(0, @(x, t, u, DuDx) diffusion_equation(x, t, u, DuDx, D1, D2, a), ...
       @diffusion_ini, @diffusion_boun, x, t);

   % Extract the solution
   u = sol(:,:,1);

   %plot the results
   figure;
   surf(x,t, u);
   xlabel('spatial coordinate (x)');
   ylabel('Time(t)');
   zlabel('Transport property profile (u)');
   title('1D diffusion model of the T cell density');

%% Provide inputs                                      vv  vv  v
function [c, f, s] = diffusion_equation(x, t, u, DuDx, D1, D2, a)

    D = zeros(size(x));  % Initialize D as a vector of zeros
    D(x <= 2/a) = D1;     % Set D1 where the condition is true
    D(x > 2/a) = D2;      % Set D2 where the condition is false

    c = 1;
    f = D .* DuDx;
    s = 0;
end
%define the initial condition
function u0 = diffusion_ini(x)
    u0 = 10^3 * ones(size(x));
end
%define the boundary condition
function [pl,ql,pr,qr] = diffusion_boun(xl,ul,xr,ur,t)
    pl = ul;
    ql = 0;
    pr = ur;
    qr = 0;
end

[채현]

I must have added the parameters D1,D2, and a when defining the diffusion_equation function. Now the graph is plotted. The problem is solved. I will now try to plot several graphs while adjusting the values of the coefficients.

Thanks for the detailed explanation!

[Dyuman Joshi]

You're welcome!