计算科学 - 数值计算平流方程 - 吾爱随笔录

数值计算平流方程

计算科学 matlab 傅立叶分析平流谱法傅里叶变换

2021-12-26 16:49:27

我正在尝试编写一个程序来计算平流方程。

u_{t} + u_{x} = 0

$u_t +u_x = 0$

我对空间导数使用谱法 $u_x$ 以及时间导数的蛙跳法 $u_t$ .

这将给出下一个时间步长的方程：

u_{j}^{k + 1} = u_{j}^{k - 1} - \frac{Δ t}{n} [f^{- 1} (i N π f (U_{j}^{k}))]

$u_{j}^{k+1} = u_{j}^{k-1} - \frac{\Delta t}{n}[f^{-1}(iN \pi f(U_{j}^{k}))]$

在哪里

$N = -n,...,n-1$ , $f =$ 傅里叶变换， $f^{-1} =$ 傅里叶逆变换， $n$ =离散点的数量。

因为我需要两个初始值（向量），所以我使用 ode45 来计算解决方案 $u_{j}^{k+1}$ 有时 $\Delta t$ 并且已经有了解决方案 $u_{j}^{k}$ 在 $t = 0$ 我们可以使用这两个向量来启动算法。

我得到错误的结果，但我看不出我做错了什么。这是一个情节

也许有人可以看到我做错了什么。

代码：

w = 0.2;
A = 1;
n = 2^8;
deltat = 2*A/n;
h = 2*A/n;   %step size for N
x = -A:h:A-h; %Discretization of N
tm = 0:deltat:2;
u = @(x) exp(-(x/w).^2);
N = -n/2:n/2-1;   

sol = zeros(length(tm),length(x)); % rows are space / columns are time;

u0 = u(x); % initial condition at t = 0
[t, ksol] = ode45(@f,[0 deltat], u0);
u1 = ksol(end,:);
sol(2,:) = u1;
sol(1,:) = u0;

for k = 3:length(tm)+1
    f_hat = fft(sol(k-1,:),n);    % fft transform
    umiddle = deltat/n * ifft(1i * pi * N * f_hat); %% inverse fft
    sol(k,:) = sol(k-2,:) - umiddle;

end

mesh(abs(real(sol(2:end,:))))

功能：

  function sol = f(~,u)
    n = 2^8;
    N = -n/2:n/2-1;   
    D = diag(1i*N*pi); % diagonal matrix (* each u0 element)
    sol = -ifft(D*fft(u,n));
    end

1个回答

我看到几个问题：

计算的 DFTfft将零模式放在数组的开头，如果要计算导数，则需要将fftshift/ifftshift应用于数组N以确保导数正确。通过用笔和纸计算出正确的表达方式，您自己很容易看到正确的表达方式，另请参阅fftshift.
您将导数项除以 $n$ ，但我不明白为什么这是正确的：解决方案中的最高频率模式在大小的网格上 $N$ 将会 $u(x) = \cos(\pi x N/(2A))$ , 有阶导数 $N$ . 然而，你计算导数的方式，最高频率模式将有一个导数 $\leq 1$ ，这似乎不对。
我认为您在导数项中也损失了两倍，因为近似于 $u_t$ 应该 $(u^{k+1}-u^{k-1})/(2\delta t)$ .
您的 $\delta t$ 太大。如果我们在冯诺依曼意义上分析您的方案的稳定性，则替换本征函数 $u(t,x) = e^{\mathrm{i}(\lambda t-\mu x)}$ ，色散关系为 $\mu\,\delta t = \sin(\lambda\,\delta t)$ . 网格上的最高频率模式 $[-A,A]$ 的 $N$ 点是 $(-1)^j = \cos(\pi N\frac{(x/A)+1}{2})$ ，所以 $\delta t$ 最多应该是 $2A/(\pi N)$ 为了稳定。

这是变化，这给了我以速度 1 传播的波的预期结果。

w = 0.2;
A = 1;
n = 2^8;
deltat = 2*A/(pi*1.001*n); %%% Try changing 1.001 to 0.999
h = 2*A/n;   %step size for N
x = -A:h:A-h; %Discretization of N
tm = 0:deltat:2;
u = @(x) exp(-(x/w).^2);
N = -n/2:n/2-1;

sol = zeros(length(tm),length(x)); % rows are space / columns are time;

u0 = u(x); % initial condition at t = 0
%[t, ksol] = ode45(@f,[0 deltat], u0);
%u1 = ksol(end,:);
sol(1,:) = u0;
u1 = u0 + deltat * ifft(-1i*pi*ifftshift(N).*fft(sol(1,:)));
sol(2,:) = u1;

for k = 3:length(tm)+1
    f_hat = fft(sol(k-1,:));    % fft transform
    umiddle = ifft(-1i * pi * ifftshift(N) .* f_hat); %% inverse fft
    sol(k,:) = sol(k-2,:) + 2*deltat * umiddle;
    u1  = sol(k-1,:);
end

[meshX,meshT] = meshgrid(x, tm);
figure;
pcolor(meshT, meshX, real(sol(2:end,:)));
shading flat;

其它你可能感兴趣的问题

上一篇这是什么PRNG功能？下一篇2d Euler 制造的解决方案