Trying to solve PDE using Finite Differential Method
2 views (last 30 days)
Show older comments
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');
0 Comments
Answers (1)
Torsten
on 15 Dec 2023
Moved: Torsten
on 15 Dec 2023
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.
4 Comments
Torsten
on 16 Dec 2023
Edited: Torsten
on 16 Dec 2023
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 ?
See Also
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!