使用PDE求解器求解PDE的高初始值u

我想使用MATLAB中的PDE工具箱解决一个简单的热方程，使用以下代码：

clear model;
width = 1e-5;
height = 1e-5;
gdm = [3 4 0 width width 0 0 0 height height]';
g = decsg(gdm, 'S1', ('S1')');
numberOfPDE = 1;
model = createpde(numberOfPDE);
geometryFromEdges(model, g);

% figure;
% pdegplot(model,'EdgeLabels','on');
% axis([-.1*width 1.1*height -.1*width 1.1*height]);
% title 'Geometry With Edge Labels Displayed';

hmax = .1 * width;
msh = generateMesh(model, 'Hmax', hmax);
p = msh.Nodes;

% figure;
% pdeplot(model);
% axis equal
% title 'Plate With Triangular Element Mesh'
% xlabel 'X-coordinate, meters'
% ylabel 'Y-coordinate, meters'

pulseDuration = 1e-13;
startTime = -2 * pulseDuration;
endTime = 5 * pulseDuration;
timeStep = 1e-2 * pulseDuration;
tlist = startTime:timeStep:endTime;    
factor_val = 1e-3;
%tlist = -1e-8 * factor_val:1e-10 * factor_val:1e-8 * factor_val;    
pulse_intensity = create_gaussian_pulse(startTime, endTime, pulseDuration, size(tlist, 2), 1e18);
figure;
plot(tlist, pulse_intensity);

f = @(region, state) (region.x < width * 1e-2) * pulse_intensity(state.t / pulseDuration);    
specifyCoefficients(model,'m',0,'d',1,'c',0.0018,'a',0,'f',@heat_source);
setInitialConditions(model, 1e18);    
numNodes = size(p, 2);
u0(1:numNodes) = 1e18;    
model.SolverOptions.RelativeTolerance = 1e-4;
model.SolverOptions.AbsoluteTolerance = 1e-5;    
R = solvepde(model, tlist);
u = R.NodalSolution;

% figure;
% plot(tlist,u(3, :));
% grid on
% title 'Temperature Along the Top Edge of the Plate as a Function of Time'
% xlabel 'Time, seconds'
% ylabel 'Temperature, degrees-Kelvin'

figure;
pdeplot(model,'XYData',u(:,end),'ColorMap','hot');
title(sprintf('Temperature In The Plate, Transient Solution( %d seconds)\n', ...
  tlist(1,end)));
xlabel 'X-coordinate, meters'
ylabel 'Y-coordinate, meters'
axis equal;

函数 heat_source 是：

function source_term = heat_source(region, state)
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here
    width = 1e-5;
%     pulseDuration = 1e-13;
%     startTime = -2 * pulseDuration;
%     endTime = 5 * pulseDuration;
%     timeStep = 1e-17;
%     tlist = startTime:timeStep:endTime;
%     pulse_intensity = create_gaussian_pulse(startTime, endTime, pulseDuration, size(tlist, 2), 1e34);
%     if(region.x < width * 1e-3)
%         source_term = 1;
%     else
%         source_term = 0;
%     end
    state.time
    source_term = (region.x < width * 1e-3) * 1;%pulse_intensity((state.time - startTime)/timeStep + 1);
end

现在我做了两个观察：

• 由于 u （ 1e18 ）的初始值较高，之后的图失败了：

错误设置类'Axes'的属性'CLim'：值必须是数字类型的1x2向量，其中第二个元素大于第一个元素并且可能是Inf

由于功能：

caxis([cmin cmax]);

其中 cmin = 1e18 和 cmax = 1e18 + 1 ，大约等于 cmin ，因此无法满足 cmin < cmax 的要求

• 我从加热函数得到的时间是 -2e-13 ， 0 或 NaN ，但不会随着时间的推移而增加，正如我所料：

ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

  -2.0000e-13


ans =

  -2.0000e-13


ans =

  -2.0000e-13


ans =

  -2.0000e-13


ans =

  -2.0000e-13


ans =

   NaN


ans =

   NaN


ans =

   NaN


ans =

   NaN


ans =

   NaN


ans =

  -2.0000e-13


ans =

  -2.0000e-13


ans =

  -2.0000e-13


ans =

  -2.0000e-13


ans =

  -2.0000e-13


ans =

     0

因此：为什么我的脚本表现如此，我该如何解决？