Using the adaptmesh Function for a Time-Dependent Source Function
显示 更早的评论
Hello !
I am trying to solve a time-dependent problem using the adaptmesh function in MATLAB, where the source term in Poisson's equation changes with time. Specifically, I want to update the source term at each time step and adapt the mesh dynamically to improve the solution’s accuracy at each stage.
The problem involves solving the Poisson equation with a time-dependent source function that has the following form:
where
is a sinusoidal function that represents the moving source:
is a sinusoidal function that represents the moving source:
Here’s the basic structure of my code, but it’s not working as expected:
% Constants
R0 = 1e-15; % Approximation of the distribution radius
L = 1e-13; % Domain size (m)
v = 1e6; % Proton speed (m/s)
tmax = 2*L/v; % Maximum time
dt = 1e-20;
A = 1e-14; % Amplitude (m)
f = 1e9; % Frequency (Hz)
phi = 0; % Initial phase (radians)
x0 = @(t) A * sin(2 * pi * f * t + phi);
model = createpde;
% Geometry (rectangle)
R1 = [3,4,-L/2,L/2,L/2,-L/2,L/4,L/4,-L/3,-L/3]'; % Rectangle
g = decsg(R1);
geometryFromEdges(model,g);
figure;
pdegplot(model,"EdgeLabels","on");
axis equal
title("Geometry With Edge Labels Displayed")
% Coefficients for Poisson's equation
c = 1; % Diffusion coefficient
a = 0; % Reaction term
sigma = 1e-15; % Standard deviation
% Time loop
t = 0;
while t <= tmax
% Update the source term at each time step (f needs to be updated)
f = @(location, state) exp(-(location.x - x0(t)).^2 / (2 * sigma^2)) / (sqrt(2 * pi) * sigma);
% Specify coefficients for Poisson's equation
specifyCoefficients(model,"m",0,"d",0,"c",c,"a",a,"f",f);
% Boundary conditions
applyBoundaryCondition(model,"dirichlet","Edge",(1:4),"u",0);
% Adapt the mesh
[u,p,e,t] = adaptmesh(g,model,c,a,f, 0, "tripick", ...
"circlepick", ...
"maxt",2000, ...
"par",1e-3);
% Solve the PDE
results = solvepde(model);
% Extract the results
u = results.NodalSolution;
% Display the results
figure;
pdeplot(model,"XYData",u, "ZData",u);
title(['Solution at time t = ', num2str(t)]);
% Update time
t = t + dt;
end
The issue is that the adaptmesh function does not accept my function f as it is written in the code. How can we dynamically adapt the mesh over time, specifically near a given location (such as the proton's position)?
Thank you in adavance for your help !
12 个评论
Why do you set the first and second order time derivatives to 0
specifyCoefficients(model,"m",0,"d",0,"c",c,"a",a,"f",f);
, but specify a time-dependent source term
f = @(location, state) exp(-(location.x - x0(t)).^2 / (2 * sigma^2)) / (sqrt(2 * pi) * sigma);
?
What is the problem you are trying to solve ? Do you want to get u == 0 as a solution for all values for t ?
Safiya
2024-11-24
Here is the initial equation I want to solve : 
And I used a gaussian function to approximate the Dirac function.
You can try to solve this equation for fixed values of t, but I'm quite sure that there should be a term in the equation that "connects" these stationary solutions over time.
Look up "Fundamental Solution" under
Here, the stationary solution for a fixed value of t is explicitly given.
Safiya
2024-11-24
@Safiya tries to solve the problem
-(d^2u/dx^2 + d^2u/dy^2) = dirac([x y]-[x0 y0])
by approximating
dirac([x y]-[x0 y0])
by
% Update the source term at each time step (f needs to be updated)
f = @(location, state) exp(-(location.x - x0(t)).^2 / (2 * sigma^2)) / (sqrt(2 * pi) * sigma);
(for whatever reason, the y-dimension is ignored)
and setting boundary conditions as u == 0 on the edges of a rectangle.
I mistakenly wrote u == 0 is the solution, but it should be
u(x,y) = -log(sqrt((x-x0)^2+(y-y0)^2))/(2*pi)
according to
(at least if the problem is defined without boundary condition).
Paul
2024-11-25
Just to make sure I'm following, the phi in the equation shown in this comment is u in the solution. But now I'm confused about
d^2u/dx^2 + d^2u/dy^2 = dirac([x y]-[x0 y0])
The left hand side is a scalar, but the right hand side is vector?
Torsten
2024-11-25
Mathematically, dirac(x) = 0 for x ~= 0 and dirac(x) = Inf for x = 0 where x can be any vector in IR^n.
Maybe "dirac" in MATLAB is differently defined.
You used this example from the UserGuide to start with ?
c = -1;
a = 0;
f = @circlef;
numberOfPDE = 1;
model = createpde(numberOfPDE);
%Create a geometry and include it in the model.
g = @circleg;
%L = 1e-13;
%R1 = [3,4,-L/2,L/2,L/2,-L/2,L/4,L/4,-L/3,-L/3]'; % Rectangle
%g = decsg(R1);
geometryFromEdges(model,g);
%Plot the geometry and display the edge labels.
figure;
pdegplot(model,"EdgeLabels","on");
axis equal
title("Geometry With Edge Labels Displayed")
applyBoundaryCondition(model,"dirichlet","Edge",(1:4),"u",0);
[u,p,e,t] = adaptmesh(g,model,c,a,f,"tripick", ...
"circlepick", ...
"maxt",2000, ...
"par",1e-3);
figure;
pdemesh(p,e,t);
axis equal
%Plot the error values.
x = p(1,:)';
y = p(2,:)';
r = sqrt(x.^2+y.^2);
uu = log(r)/2/pi;
figure;
pdeplot(p,e,t,"XYData",u-uu,"ZData",u-uu,"Mesh","off");
figure;
pdeplot(p,e,t,"XYData",u,"ZData",u,"Mesh","off");
Safiya
2024-11-25
Yes, that's right. I tried to adapt this code to my case where the refinement position of the mesh is updated with the position of the proton moving in time (
from 
).
from ).For your code so far, you can just plot the result as
pdeplot([p(1,:)+x(time);p(2,:)+y(time)],e,t,"XYData",u-uu,"ZData",u-uu,"Mesh","off");
where [x(time),y(time)] is the position of the proton at t = time.
This works since the results for different values of t don't depend on one another and because the solution for (0,0) is just translated by the vector (x(time),y(time)).
回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Geometry and Mesh 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
