关于在MATLAB中模拟抛物型PDE解的爆破现象及温和解的技术问询

Great questions! Let's break down how to approach both your simulation goals in MATLAB, along with practical tips and learning resources.

一、模拟Fujita方程的爆破/全局解行为 #

You're right that we can't simulate the infinite domain $\mathbb{R}^n$ directly, but we can use large bounded domains with appropriate boundary conditions to approximate the infinite-domain behavior. Here's how to do it:

1. Using MATLAB's PDE Toolbox #

The PDE Toolbox is a beginner-friendly starting point for 1D/2D cases:

• Define a large rectangular or circular domain (large enough that the solution hasn't reached the boundary by the time you expect blow-up or steady-state).
• Set boundary conditions that mimic the infinite domain: for example, use Dirichlet boundary conditions set to 0 (since typical initial conditions like Gaussian functions decay to 0 at infinity).
• Define the parabolic PDE and run time-dependent simulations:

  % Example for 2D; adjust for 1D as needed
  pde = createpde();
  geometryFromEdges(pde,@squareg); % Use a large square, e.g., [-10,10]x[-10,10]
  p = 1.5; % Test blow-up case for n=2 (critical exponent is 2)
  specifyCoefficients(pde,'m',0,'d',1,'c',1,'a',0,'f',@(region,state) state.u.^p);
  setInitialConditions(pde,@(location) exp(-(location.x.^2 + location.y.^2))); % Gaussian initial condition
  tlist = linspace(0,5,100);
  results = solvepde(pde,tlist);
  

• Monitor the maximum value of $u$ to detect blow-up:

  max_u = arrayfun(@(k) max(results.NodalSolution(:,k)),1:length(tlist));
  plot(tlist,max_u,'LineWidth',1.2);
  xlabel('Time'); ylabel('Max u(t)');
  

  If max_u spikes to a very large value (or causes numerical overflow), that's your numerical indication of finite-time blow-up.

2. Custom Finite Difference/Finite Element Code #

For more control (especially for 3D or custom boundary conditions), write your own code:

• Use a uniform spatial grid and implicit time-stepping (like Crank-Nicolson) to avoid stability issues with explicit methods.
• For the Laplacian, use second-order central differences in 1D/2D.
• Iterate over time steps, updating $u$ each iteration, and track $|u(t)|_\infty = \max(u(:))$. When this value exceeds a pre-defined threshold (e.g., $10^6$), you can conclude blow-up has occurred numerically.

Key Note: Test with $p$ values on either side of the critical exponent $1+2/n$ (e.g., $n=2$, $p=1.5$ for blow-up, $p=2.5$ for global solutions) to verify the behavior matches theory.

二、模拟带$|x|^\gamma u^p$的温和解 #

Mild solutions are defined via integral equations, so we can use time-stepped convolution methods with the heat semigroup $S(t)$. Here's a practical approach:

1. Core Idea: Discretize the Integral Equation #

We can approximate the integral in the mild solution using a time-step sum:
$$u(t_n) \approx S(t_n)u_0 + \sum_{k=0}^{n-1} S(t_n - t_k) \left(|x|^\gamma u^p(t_k)\right) \Delta t$$
The critical part is efficiently computing $S(t)f$, which is the convolution of $f$ with the heat kernel $K(x,t)$.

2. Accelerate Convolutions with FFT #

Direct convolution is slow for large grids, so use Fast Fourier Transforms (FFT) to compute $S(t)f$ in the frequency domain:

• Compute the FFT of the heat kernel $K(x,t)$ for each time step.
• Compute the FFT of your function $f$, multiply by the FFT of $K(x,t)$, then take the inverse FFT to get $S(t)f$.

3. MATLAB Implementation Steps #

• Define a large spatial grid $x$ (e.g., $x \in [-10,10]^n$) and precompute $|x|^\gamma$ (note: for $\gamma<0$, add a small epsilon like $1e-6$ to $|x|$ to avoid division by zero at $x=0$).
• Initialize $u_0$ (e.g., a Gaussian function).
• Precompute the heat kernel's FFT for each time interval $t_n - t_k$.
• Iterate over time steps:
  1. Compute $S(t_n)u_0$ via FFT convolution.
  2. For each previous time step, compute $|x|^\gamma \cdot u(t_k)^p$, apply $S(t_n-t_k)$ via FFT, multiply by $\Delta t$, and accumulate the sum.
  3. Update $u(t_n)$ and monitor its maximum value.

4. Handling Singularities #

For $\gamma<0$, $|x|^\gamma$ diverges at $x=0$. You can:

• Use a non-uniform grid that's denser near $x=0$ to capture the singularity.
• Regularize the singularity by replacing $|x|$ with $\sqrt{|x|^2 + \epsilon^2}$ for small $\epsilon>0$.

三、推荐学习书籍 #

• Numerical Solutions of Partial Differential Equations by K.W. Morton and D.F. Mayers: A classic, accessible introduction to finite difference/finite element methods for PDEs, including parabolic equations.
• Partial Differential Equations: An Introduction by Walter A. Strauss: Balances theory and intuition, with examples of blow-up phenomena and numerical approaches.
• Numerical Methods for Nonlinear Partial Differential Equations by J.W. Thomas: Dives deeper into nonlinear PDEs, including techniques for handling integral formulations like mild solutions.

备注：内容来源于stack exchange，提问作者Ilovemath