How to solve the following PDE equation

Answer 1

[Torsten]

I think you mean x(s,0) = 10 instead of x(s,t) = 10, don't you ?

The easiest way to solve the equation is to discretize the expressions on the right-hand side and solve the resulting system of ordinary differential equations using ode15. Here, the integral term can be approximated by MATLAB's "trapz".

Look up "method-of-lines" for more details.

[Sam]

I tried the solving it using pdepe, but was getting this error

"Error using matlab.internal.math.getdimarg

Dimension argument must be a positive integer scalar within indexing range."

The code is attached below

function theta = solvePDE

x = linspace(0, 1, 50);

t = linspace(0, 10, 100);

sol = pdepe(0, @pdefun, @icfun, @bcfun, x, t);

theta = sol(:,:,1);

end

function [c, f, s] = pdefun(x, t, u, dudx)

integral_term = trapz(x, u);

c = 1;

f = dudx;

s = integral_term;

end

function u0 = icfun(x)

u0 = 30* ones(size(x));

end

function [pl, ql, pr, qr] = bcfun(xl, ul, xr, ur, t)

pl = ul - 30;

ql = 0;

pr = ur - 30;

qr = 0;

end

[Torsten]

integral_term = trapz(x, u);

"pdepe" supplies x and u pointwise, not for the complete interval [0,1]. Thus using "pdepe" this way is not possible.

It might be possible using this code for your purpose:

https://uk.mathworks.com/matlabcentral/fileexchange/97437-pde1dm

Answer 2

[Bill Greene]

I don't know how to evaluate that integral using pdepe. However, I have written a pde solver (pde1dm) that has an input syntax very similar to pdepe and includes an option that makes this straightforward.

The "vectorized" option tells pde1dm to call your pdefun with a vector of x values spanning the complete spatial domain of your problem. I have included a slightly-modified version of your code below. If you want to try pde1dm, it can be downloaded using the link above.

function matlabAnswers_7_12_2024
theta=solvePDE;
end

function theta = solvePDE
nx=50;
x = linspace(0, 1, nx);
nx2=ceil(nx/2);
t = linspace(0, 10, 100);
if 0
    sol = pdepe(0, @pdefun, @icfun, @bcfun, x, t);
else
    opts.vectorized='on';
    sol = pde1dm(0, @pdefun, @icfun, @bcfun, x, t,opts);
end
theta = sol(:,:,1);
figure; plot(x, sol(end,:)); title 'solution at final time';
figure; plot(t, sol(:,nx2)); title 'solution at center as a function of time';
end

function [c, f, s] = pdefun(x, t, u, dudx)
nx=length(x);
integral_term = trapz(x, u);
c = ones(1,nx);
f = dudx;
s = ones(1,nx)*integral_term;
end

function u0 = icfun(x)
u0 = 30* ones(size(x));
end

function [pl, ql, pr, qr] = bcfun(xl, ul, xr, ur, t)
pl = ul - 30;
ql = 0;
pr = ur - 30;
qr = 0;
end