b.c. issue

xzczd

11:27

囧，你一次开了两个房间

——今晚应该都会在线，不过有点事要折腾，可能回复不会很及时

xiaohuamiao

11:58

Hi, sorry to bother you again.
When I use your FDM code to solve ODEs, I recently noticed that, in some cases, the specified boundary conditions cannot be reached in the end. A MWE is as follows.
I only want vanishing Dirichlet right bc. And I neglect the left bc since it is 1st order. You can see the right bc is not satisfied in general.
eps0 = 1*^-10; domain = {eps, Rcutoff}; points = 600; difforder = 4; lbcvalue = 1; rbcvalue = 0;
m = 1.0; A[l_, r_] := l/r;
Fop1[y_, l_, pm_] := I (-D[y, r] + pm A[l, r] y);

哈你好啊已经用英文打了就这么发了

@xzczd l是任意整数

@xzczd 我只解过1st order和2nd order 目前只在1st order里发现这个问题不清楚是不是1st order的某种特性

xiaohuamiao

12:18

eigenvalue problem对wrong Dirichlet bc不敏感这件事大体是确实的(我没试过非Dirichlet的) 不过我没见到什么地方提及这个试着简单证明过貌似也不是很好搞…… 但是不论如何真要是fix住了boundary value 那解在边界上貌似也不应该是这个小代码画出来的那样所以Q2想问你这个code实际是怎么fix bc的

xiaohuamiao

12:46

@xzczd 我艹抱歉了你可以先忽略这个具体的code跟Q1 这个问题有可能恰好是bc满足不了的我发现我把有个相关的数学定理给记错了…… 所以主要是问Q2

xzczd

13:13

网突然很渣……

你这个问题为什么对边界不敏感我也一直没搞清楚……这问题有没有解析解啊，有的话从那边着手分析可能会好一点。

引入边界的思路其实一直是一样的啊，就是把边界的方程给混进方程组里去

只是平时解微分方程的时候，对方程形式要求更宽松，所以一惯都是删掉边界上的方程给边界条件让出位置

但是求特征值的时候没法直接删，所以用的是这种代入的方法

xiaohuamiao

“是求特征值的时候没法直接删，所以用的是这种代入的方法”
啊？什么意思我怎么没看出来有什么不同

xzczd

13:28

因为特征值问题对方程形式有要求啊，一侧必须是 特征值 × y[x]的形式

……等一下，对于纯Dirhchlet边界或许可以换一种代入的方法。lhsL 和 lhsR 就不用定义了，直接把边界代入到ae里面去，效果（应该）是一样的。这样编程会简单一点点。@xiaohuamiao

xiaohuamiao

我一直以为边界上的值是设定好了的解的时候不会变为什么有时候解画出来边界上的值就会变的跟设定的不一样了呢

xzczd

直接删掉几个方程再拿边界塞上去的话就丧失了上面那种形式。

@xiaohuamiao 如果出现了这种现象，那应该就是我上面说的原因了。因为现在这个方法不是完全的固定，而是部分的替代。

……等下，好像不对

会不会是和EigenVector 这个函数有关系……

啊啊啊特征值的求解原理我搞不清楚啊啊啊

但是应该和解方程是不一样的，吧

xiaohuamiao

两件事
1 设定好的这个值居然会变
2 错的会自动变成对的（实际使用上经常貌似是件好事……）这肯定是特征值问题的一般特性虽然一下子并不好证明

我现在主要是想理解1……

xzczd

13:56

@xiaohuamiao 刚试了一下，发现对比较大的阵，只改遍其中一行的取值，对其特征值的绝对值的影响是可以忽略的。不知道这个是不是原因？

n = 100;
mat = RandomReal[1, {n, n}];
mat2 = Block[{mat = mat}, mat[[1]] = RandomReal[1, n]; mat];
lst = Eigenvalues /@ {mat, mat2};

ListLinePlot[lst // Abs, PlotRange -> All]
ListLinePlot[lst // Im, PlotRange -> All]

你可以试试

当然，实部和虚部分开画的话就明显了

感觉应该是部分的原因，吧……

xiaohuamiao

14:09

@xzczd 我这个小code是照你答案里那样写的我理解就是删了几个方程然后直接加进去边界条件你为什么说“求特征值的时候没法直接删，所以用的是这种代入的方法” 我哪里有误解吗

xzczd

14:31

代入发生在 e.g. lhsR = lhs/. ……这步

@xiaohuamiao 我说的“没有直接删”指的是，这里没有直接用边界条件的方程代替原方程，而是先把边界条件混入原方程以形成新方程，再把这个新方程拿去替换。

但我现在觉得这点区别应该和特征值的问题没（什么大的）关系

@xiaohuamiao 我说的“没有直接删”指的是，这里没有直接用边界条件的方程代替原方程，而是先把边界条件混入原方程以形成新方程，再把这个新方程拿去替换。

Transcript for