Trying to solve PDE using Finite Differential Method

I am trying to solve the follwing equations using the finite differential method. I got partway through but it seems to be outputting zero for x, y, and z. Can I get help on what is going wrong?

dxdt=d1*d^2x/dl^2+r*x*(1-x/K)-a*x-w*x

dydt=d2*d^2y/dl^2g*w*x-c*y-b*y*z

dzdt=d3*d^2z/dl^2f*b*y*z-s*z

% Define parameters

r = 1.3;

K = 700;

a = 0.2;

w = 0.53;

c = 0.001;

b = 0.015;

g = 1;

s = 0.7;

f = 1;

d1 = 0.01;

d2 = 0.03;

d3 = 0.000;

% Define grid and step sizes

N = 100; % Number of grid points

L = 10; % Length of the domain

dl = L / N;

l = linspace(0, L, N);

% Initial conditions

x0 = zeros(1, N);

y0 = zeros(1, N);

z0 = zeros(1, N);

% Solve the PDEs using finite differences

for t = 1:100 % Time steps

% Calculate second derivatives using finite differences

d2xdl2 = (circshift(x0, -1) - 2 * x0 + circshift(x0, 1)) / dl^2;

d2ydl2 = (circshift(y0, -1) - 2 * y0 + circshift(y0, 1)) / dl^2;

d2zdl2 = (circshift(z0, -1) - 2 * z0 + circshift(z0, 1)) / dl^2;

% Updating x, y, z using the PDEs

dxdt = d1 * d2xdl2 + r * x0 .* (1 - x0 / K) - a * x0 - w * x0;

dydt = d2 * d2ydl2 + g * w * x0 - c * y0 - b * y0 .* z0;

dzdt = d3 * d2zdl2 + f * b * y0 .* z0 - s * z0;

% Time integration using forward Euler's method

x = x0 + dl^2 * dxdt;

y = y0 + dl^2 * dydt;

z = z0 + dl^2 * dzdt;

% Update for the next time step

x0 = x;

y0 = y;

z0 = z;

end

% Plot the solutions

figure;

subplot(3, 1, 1);

plot(l, x);

xlabel('l');

ylabel('x');

title('Solution for x');

subplot(3, 1, 2);

plot(l, y);

xlabel('l');

ylabel('y');

title('Solution for y');

subplot(3, 1, 3);

plot(l, z);

xlabel('l');

ylabel('z');

title('Solution for z');

Answer 1

[Torsten]

What are your boundary conditions for x, y and z at l = 0 and l = L ?

In your code, you seem to use initial conditions and boundary conditions equal to 0 for x, y and z. So the solution x = y = z = 0 for all l and all t is correct.

[Daniel Mai]

Thanks. I forgot to add conditions. I have modified it slightly to add the conditions and to simplify my equations.

solvePDE()

function solvePDE()
% Define parameters
r = 1.3;
K = 700;
a = 0.2;
w = 0.53;
c = 0.001;
b = 0.015;
g = 1;
s = 0.7;
f = 1;
d1 = 0.01;
d2 = 0.03;
d3 = 0.015;

% Define spatial grid
N = 110; % Number of spatial points (including extra points for boundary)
L = 110; % Length of the domain (including extra length for boundary)
dl = L / (N - 1);
l = linspace(-10, 100, N);

% Initial conditions
x = 1000 * ones(1, N); % Initial condition for x (baby turtle population)
y = 100 * ones(1, N);  % Initial condition for y (juvenile turtle population)
z = ones(1, N);        % Initial condition for z (predator population)

% Apply boundary conditions
z(end) = 200; % Boundary condition at l = 100
z(1) = 10;    % Boundary condition at l = -10

% Update initial condition for z (predator population)
z = updateInitialConditionZ(z, N);

% Time parameters
timesteps = 100; % Number of time steps
dt = 0.1;        % Time step size

% Implement finite difference method
for t = 1:timesteps
    % Calculate second derivatives using central differences for dz/dl^2
    d2zdl2 = (circshift(z, -1) - 2 * z + circshift(z, 1)) / (dl^2);
    
    % Update z using the PDE and boundary conditions
    dzdt = d3 * d2zdl2 .* (f * b * y .* z - s) - s * z;
    z = z + dt * dzdt;
    
    % Update x using the differential equation dx/dt
    dxdt = r * x .* (1 - x / K) - a * x - w * x;
    x = x + dt * dxdt;
    
    % Update y using the differential equation dy/dt
    dydt = g * w * x - c * y - b * y .* z;
    y = y + dt * dydt;
end

% Plot the solutions
figure;
subplot(3, 1, 1);
plot(l, x);
xlabel('Distance from shore (l)');
ylabel('Baby Turtle Population (x)');
title('Population Dynamics: Baby Turtles');

subplot(3, 1, 2);
plot(l, y);
xlabel('Distance from shore (l)');
ylabel('Juvenile Turtle Population (y)');
title('Population Dynamics: Juvenile Turtles');

subplot(3, 1, 3);
plot(l, z);
xlabel('Distance from shore (l)');
ylabel('Predator Population (z)');
title('Population Dynamics: Predators');
end

function z = updateInitialConditionZ(z, N)
% Update initial condition for z (predator population)
initialPredator = 20; % Initial predator population at shore
maxPredator = 200; % Maximum predator population
gradient = (maxPredator - initialPredator) / (N - 11); % Gradient

for i = 1:N
    if i >= 2 && i <= 101
        z(i) = initialPredator + (i - 1) * gradient;
    else
        z(i) = 10; % Setting constant value for boundary regions
    end
end
end

[Torsten]

So problem solved ?

[Daniel Mai]

Think so?

[Torsten]

You can check using MATLAB's "pdepe".

I think "circshift" is only the correct tool to build the second-order derivatives from the z-vector if the boundary condition is a periodic one. Is this the case for your problem ?

And in the code above, you neglect diffusion for x and y ?