标签: pde

为什么(A + B)的FFT与FFT(A)+ FFT(B)不同?

我差不多一个月一直在和一个非常奇怪的虫子打架.问你们这是我最后的希望.我在C中编写了一个程序,它使用傅里叶(或倒数)空间中的隐式欧拉(IE)方案集成了2d Cahn-Hilliard方程:

IE方法

"帽子"意味着我们处于傅里叶空间:h_q(t_n + 1)和h_q(t_n)是时间t_n和t_(n + 1)的h(x,y)的FT,N [h_q]是非线性算子应用于傅立叶空间中的h_q,而L_q是线性的,同样在傅立叶空间中.我不想过多介绍我使用的数值方法的细节,因为我确信问题不是来自那里(我尝试使用其他方案).

我的代码实际上非常简单.这是开始,基本上我声明变量,分配内存并为FFTW例程创建计划.

# include <stdlib.h>
# include <stdio.h>
# include <time.h>
# include <math.h>
# include <fftw3.h>
# define pi M_PI

int main(){

// define lattice size and spacing
int Nx = 150;         // n of points on x
int Ny = 150;         // n of points on y
double dx = 0.5;      // bin size on x and y

// define simulation time and time step
long int Nt = 1000; …
Run Code Online (Sandbox Code Playgroud)

c fftw pde dft

49
推荐指数
1
解决办法
898
查看次数

Eclipse PDT中的"API基线"是什么

自升级到Eclipse 3.7以来,Eclipse PDE插件希望我为所有Eclipse插件项目指定"API Baseline".

然而,似乎没有文件实际上解释了什么'API Baseline'在这里代表什么,以及它用于什么.

有人可以解释一下吗?

java eclipse eclipse-pde eclipse-plugin pde

15
推荐指数
2
解决办法
4326
查看次数

在Python中使用隐式Euler解决PDE-错误的输出

我将尝试解释发生的情况和问题。

这有点数学,因此不支持乳胶,因此遗憾的是我不得不求助于图像。我希望可以。

enter image description here

enter image description here

我不知道为什么要倒转,对此感到抱歉。无论如何,这是一个线性系统Ax = b,我们知道A和b,因此我们可以找到x,这是我们在下一个时间步长上的近似值。我们将继续执行直到时间t_final。

这是代码

import numpy as np

tau = 2 * np.pi
tau2 = tau * tau
i = complex(0,1)

def solution_f(t, x):
    return 0.5 * (np.exp(-tau * i * x) * np.exp((2 - tau2) * i * t) + np.exp(tau * i * x) * np.exp((tau2 + 4) * i * t))

def solution_g(t, x):
    return 0.5 * (np.exp(-tau * i * x) * np.exp((2 - tau2) * i * t) - np.exp(tau * i * …
Run Code Online (Sandbox Code Playgroud)

python algorithm numerical-methods pde

14
推荐指数
2
解决办法
293
查看次数

如何在R中使用外推栅格

我试图缩减使用的方法在气候条件文章使用R软件.我几乎在那里,但我错过了几个步骤

需要的包和数据

在本例中,我将一些数据上传到archive.org网站,以加载本例中使用的所需包和数据,使用以下代码:

library(raster)
library(rgdal)

download.file("https://archive.org/download/Downscaling/BatPatagonia.rds", "Bat.rds")
download.file("https://archive.org/download/Downscaling/TempMinPatNow.rds", "Tmin.rds")

BatPatagonia <- readRDS("Bat.rds")
TempMinPatNow <- readRDS("Tmin.rds")
Run Code Online (Sandbox Code Playgroud)

BatPatagonia是一个栅格文件,其中Bathymetry和从GEBCO数据集中提取和转换的区域的高度,而TempMinPatNow是从worldclim中提取的1月相同区域的最低温度.下面是数据集的图:

在此输入图像描述

这个问题的目标

为了从过去的冰川最大值下调过去的数据,我需要模拟当前的气候如果海平面与过去相同的情况.为了做到这一点,我使用GEBCO数据,并且或多或少地弄清楚海岸是什么.根据上面引用的文章中的方法,这是前面的三个步骤:

  1. 为海拔20米以上的土地创建DEM
  2. 在移动窗口中计算多元线性回归
  3. 将系数外推到海洋

第3点是我一直在努力做的事情,我将展示我如何做到前2点,并展示我一直在寻找的尝试解决第3点的问题

1.为海拔20米的土地建立DEM

为了做到这一点,我采用了BatPatagonia栅格,并使用以下代码将所有超过20米的值替换为NA值:

Elev20 <- BatPatagonia

values(Elev20) <- ifelse(values(Elev20) <= 20, NA, values(Elev20))
Run Code Online (Sandbox Code Playgroud)

生成的栅格如下图所示

在此输入图像描述

2.在移动窗口中计算多元线性回归

根据第2591页的手稿,下一步是在移动窗口中使用以下公式对超过20米的高度进行多元线性回归:

在此输入图像描述

我们已经有高程数据,但我们还需要纬度和经度的栅格,为此我们使用以下代码,首先创建纬度和经度栅格:

Latitud <- BatPatagonia
Longitud <- BatPatagonia

data_matrix <- raster::xyFromCell(BatPatagonia, 1:ncell(BatPatagonia))

values(Latitud) <- data_matrix[, 2]
values(Longitud) <- data_matrix[, 1]
Run Code Online (Sandbox Code Playgroud)

我们将乘以高度超过20米的区域的光栅掩模,这样我们只得到我们需要的值:

Elev20Mask <- BatPatagonia

values(Elev20Mask) <- ifelse(values(Elev20Mask) <= 20, NA, 1)

Longitud <- Elev20Mask*Longitud

Latitud <- Elev20Mask*Latitud
Run Code Online (Sandbox Code Playgroud)

现在我将使用响应变量和预测变量构建一个堆栈:

Preds <- …
Run Code Online (Sandbox Code Playgroud)

r pde r-raster

12
推荐指数
1
解决办法
397
查看次数

数值不稳定FFTW <> Matlab

我试图用伪谱方案数值求解Swift-Hohenberg方程http://en.wikipedia.org/wiki/Swift%E2%80%93Hohenberg_equation,其中线性项在傅立叶空间中隐式处理,而在现实空间中评估非线性.简单的Euler方案用于时间积分.
我的问题是我提出的Matlab代码完美无缺,而依赖FFTW进行傅里叶变换的C++代码变得不稳定,经过几千个步骤后就会发散.我已经跟踪了非线性项的处理方式(参见C++代码中的注释).如果我只使用Phi的真实部分,就会发生不稳定.然而,由于数值舍入误差,Phi应该只有一个可忽略的虚部,并且Matlab正在做类似的事情,保持Phi纯粹真实.在Octave下,Matlab代码也运行良好.初始条件可能类似于
R=0.02*(rand(256,256)-0.5);
Matlab(小振幅波动).

TLDR;

为什么这些代码片段做了不同的事情?具体来说,我如何使C++代码的工作方式与Matlab版本相同?

编辑1:

为了完整起见,我使用FFTW提供的R2C/C2R功能添加了代码.有关详细信息,请参阅http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html(我希望我的数据布局正确).此代码始终显示约3100个时间步后的不稳定性.如果我将dt减少到例如0.01,则会发生10次.

使用复杂DFT的C++代码

#include <iostream>
#include <fstream>
#include <cmath>
#include <fftw3.h>

int main() {

const int N=256, nSteps=10000;
const double k=2.0*M_PI/N, dt=0.1, eps=0.25;

double *Buf=(double*)fftw_malloc(N*N*sizeof(double));
double *D0=(double*)fftw_malloc(N*N*sizeof(double));

// complex arrays
fftw_complex *Phi=(fftw_complex*)fftw_malloc(N*N*sizeof(fftw_complex));
fftw_complex *PhiF=(fftw_complex*)fftw_malloc(N*N*sizeof(fftw_complex));
fftw_complex *NPhiF=(fftw_complex*)fftw_malloc(N*N*sizeof(fftw_complex));

// plans for Fourier transforms
fftw_plan phiPlan=fftw_plan_dft_2d(N, N, Phi, PhiF, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_plan nPhiPlan=fftw_plan_dft_2d(N, N, NPhiF, NPhiF, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_plan phiInvPlan=fftw_plan_dft_2d(N, N, Phi, Phi, FFTW_BACKWARD, FFTW_ESTIMATE);

std::ifstream fin("R.dat", std::ios::in | std::ios::binary); // read initial …
Run Code Online (Sandbox Code Playgroud)

c++ matlab numerical fft pde

9
推荐指数
1
解决办法
1463
查看次数

Eclipse PDE:什么是"目标平台"?

当您使用Eclipse PDE(插件开发环境)时,有一个术语"目标平台".这究竟是什么意思?

eclipse pde target-platform

7
推荐指数
1
解决办法
3813
查看次数

PDE的"更新站点向导"仍然是创建更新站点的正确方法吗?

这个问题的关键是从熟悉Eclipse安装系统当前状态的人那里获得验证.

我有一个Eclipse插件,我想使用最简单(但正确)的方法为它创建一个更新站点.我最初的印象是我:

  1. 使用PDE 功能项目创建功能.
  2. 将我的插件添加到该功能.
  3. 使用PDE 更新站点项目创建更新站点.
  4. 将我的功能添加到更新站点.

然后我开始在整个Eclipse安装系统上寻找文档,并开始阅读这个P2的东西,这显然是新的,取代之前的任何东西.

我找到的关于PDE项目和向导的信息没有讨论P2,或者我没有找到正确的信息,这使我有点担心我可能做错了什么或遗漏了一些重要的事情.

那么,我是否需要关心P2,或者P2是否可以安全地忽略它,只要我继续使用PDE工具?

谢谢!

eclipse eclipse-pde pde

7
推荐指数
1
解决办法
1147
查看次数

如何在eclipse中获取工作空间路径?

我正在研究PDE(Eclipse插件项目).

我需要获得工作空间路径.

我的文本小部件(swt)应该设置当前工作空间路径.

如何在eclipse中获取工作空间路径?

java eclipse workspace eclipse-plugin pde

7
推荐指数
2
解决办法
2万
查看次数

在具有3个边界条件的二维阵列中迭代有限差分格式

在完成本科毕业论文的中间步骤中,我试图通过显式有限差分方案来解决R中的Black-Scholes PDE.到目前为止,我一直在及时(沿时间离散化)和股票价格离散化反复迭代.我已经实现了欧洲看涨期权的边界条件,并且我试图根据Phil Goddard对该模型的描述来向后解决问题.

R <- matrix(NA, ncol=length(t), nrow=length(S))

# Value at Maturity:
R[,length(t)] <- Vold <- pmax(S - K,0) 

# Value at S = S_max; V = fwd
R[length(S),] <- V_max <- S_max - K * exp(-r_yr * (t_T - t))

# Value at S = S_0; V = 0
R[1,] = V_min <- t - t

# Run Model
# -----------------------------------------------------------------------------------------

# Define Coefficient-Functions a_j, b_j, and c_j
a_ <- function(j){
  a_j <- (1/2) * …
Run Code Online (Sandbox Code Playgroud)

arrays for-loop r numerical-methods pde

6
推荐指数
0
解决办法
1342
查看次数

3D 扩散/热方程的有限差分法

我正在尝试使用有限差分来求解 3D 扩散方程。我认为我的主循环有问题。特别地,离散方程为:

使用诺依曼边界条件(仅以一个面为例):

现在的代码:

import numpy as np
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D ##library for 3d projection plots
%matplotlib inline

kx = 15     #Number of points
ky = 15
kz = 15
largx = 90  #Domain length.
largy = 90
largz = 90   

dt4 = 1/2 #Time delta (arbitrary for the time).
dx4 = largx/(kx-1)    #Position deltas.
dy4 = largy/(ky-1)
dz4 = largz/(kz-1) 

Tin = 25    #Initial temperature
kapp = 0.23

Tamb3d = 150 #Ambient …
Run Code Online (Sandbox Code Playgroud)

python physics numpy heat pde

6
推荐指数
1
解决办法
7095
查看次数