MATLAB的答案

我怎么能让bvp4c解决方案更稳定?

29日视图(30天)
JK
JK 2018年2月6日
评论道: 的孩子叫K G2023年5月31日
你好,社区,
我有一个问题的稳定BVP4C-solution当导数变得很小。方程系统描述了四个组件的摩尔流沿着膜。每个物种的梯度沿膜的产品各自的渗透率Q和高压侧之间的分压差(滞留物)和低压侧(渗透)。出现问题时快速扩散的分压组件双方或多或少是相等的。
我解决边值问题的形式
y ' = dy / dx = f (y,哟)
摩尔流的y是一个向量的4种滞留物膜。在x = 0,我必须解决一个隐式方程系统归来的一组参数
g(哟,归来)= 0
并计算导数
y ' (x = 0) = f(哟,归来)
你是y (x = 0),一个未知的参数(我的左边界值——滞留物摩尔流末尾的膜)符合bvp4c-solver。我知道正确的边值在x = 1饲料进入膜滞留物一侧。
ODE非常简单:
函数就要= ydot (x, y,哟,n)
Q = [2.3369 e-08;1.17095 e-06;3.81087 e-08;2.36619 e-07] * n;
公关= 1 e6;
页= 1 e5;
如果x = = 0
选择= optimset (“显示”,“没有”);
归来= fsolve(@(归来)PermOmega(归来),Q /笔(Q)选项);
就要=问:* (y /笔(y) * pR-yPo *页);
其他的
就要=问:* (y /笔(y) *公关- (y-yo) / sum (y-yo) *页);
结束
函数g = PermOmega(归来)
%渗透成分x = 0
g = yPo-Q。*(哟/笔(哟)* pR-yPo *页)/笔(问:*(哟/笔(哟)* pR-yPo *页));
结束
结束
我最初想的解决方案之间的直线筛孔尺寸的两个边界值10。的问题是在第二行代码,渗透率的增加膜的数量n, n = 1的解决方案是稳定的(膜),但当我分压增加膜的数量变得相似,ODE变得不稳定(见图片)
我能做些什么来提高稳定性当梯度膜的小?
到目前为止我已经提供了雅可比矩阵(使用隐函数微分在x = 0),这使它速度变快,但不会更准确。我已经改变了相对和绝对误差BVP4C-solver的公差。我有试过BVP5C。我归结我原来ODE-system与三个阶段每个7变量(在滞留物和渗透压力加上滞留物温度)基本系统只有4变量。
我需要稳定的原因是,我想要使用的解决方案在一个FMINCON-optimizer有时会选择参数会导致不稳定的解决方案。金宝搏官方网站
什么好主意吗?

在回答这个问题。

接受的答案

Torsten
Torsten 2018年2月9日
编辑:Torsten 2018年2月9日
谢谢你的详细的解释这个问题。
我用下列问题公式化COMSOL多重物理量,它没有问题:
解决变量:
y,哟
方程:
就要=问:* (y /笔(y) *公关- (y-yo) / sum (y-yo) *页)
dyodx = 0
初始条件:
y = yo0 + x * (yalpha-yo0)
哟= yo0
边界条件:
哟(x = 0) = y
y (x = 1) = yalpha
请注意,我解出8未知函数,不仅4。
注意进一步ydot只在内部网格点,称为x = 0。因此除零不会发生。
最好的祝愿
Torsten。

更多的答案(6)

Torsten
Torsten 2018年2月7日
你必须指定滞留物摩尔流边界常规的x = 0 bvp4c,不是在你指定的常规ODE。
最好的祝愿
Torsten。
JK
JK 2018年2月7日
嗨Torsten,
谢谢你的回答。这里的代码调用BVP4C:
%初始条件在x = 1(滞留物入口)
yalpha = (0.3; 0.25; 0.01; 0.01);
%估计滞留物值的膜(x = 0)
yo0 = (0.95; 0.60; 0.92; 0.70)。* yalpha;
%边值解算器
FJac = @ (x, y,哟)Jacfun (x, y,哟,n);
BCJac ={[0眼(4);(4)],[0(4);眼(4)],[黑眼圈(4);0 (4)]};
bcfun = @(是的,yb,哟)[ya-yo;yb-yalpha];
initguess = @ (x) yo0 + x * (yalpha-yo0);
odefun = @ (x, y,哟)ydot (x, y,哟,n);
solinit = bvpinit (linspace (0, 1, 11), initguess, yo0);
选择= bvpset (“统计数据”,“上”,“RelTol”1 e - 3,“FJacobian”FJac,“BCJacobian”,BCJac);
索尔= bvp4c (odefun bcfun、solinit选项);
Jacfun计算两个雅可比矩阵 dy ' / dy 和 dy ' / dyo (代码未显示)每个矩阵必须分开计算x = 0到x > 0。
bcfun = @(是的,yb,哟)[ya-yo;yb-yalpha]计算边界条件。你是参数计算的边界值解算器。
我指定 y (x = 0) 在我和边界条件指定 dy / dx (x = 0) ODE-equation。我仍然不知道如何解决这个问题。我想这与ODE的刚度。
亲切的问候
JK
1评论
Torsten
Torsten 2018年2月8日
试一试
%初始条件在x = 1(滞留物入口)
yalpha = (0.3; 0.25; 0.01; 0.01);
%估计滞留物值的膜(x = 0)
yo0 = (0.95; 0.60; 0.92; 0.70)。* yalpha;
%边值解算器
bcfun = @公元前(是的,yb,哟)(是的,yb,哟,n, yalpha);
initguess = @ (x) yo0 + x * (yalpha-yo0);
odefun = @ (x, y,哟)ydot (x, y,哟,n);
solinit = bvpinit (linspace (0, 1, 11), initguess, yo0);
选择= bvpset (“统计数据”,“上”,“RelTol”1 e - 3);
索尔= bvp4c (odefun bcfun、solinit选项);
函数就要= ydot (x, y,哟,n)
Q = [2.3369 e-08;1.17095 e-06;3.81087 e-08;2.36619 e-07] * n;
公关= 1 e6;
页= 1 e5;
就要=问:* (y /笔(y) *公关- (y-yo) / sum (y-yo) *页);
结束
函数res = bc(是的,yb,哟,n, yalpha)
Q = [2.3369 e-08;1.17095 e-06;3.81087 e-08;2.36619 e-07] * n;
公关= 1 e6;
页= 1 e5;
选择= optimset (“显示”,“没有”);
归来= fsolve(@(归来)PermOmega(归来),Q /笔(Q)选项);
res = [ya-yPo; yb-yalpha];
函数g = PermOmega(归来)
%渗透成分x = 0
g = yPo-Q。*(哟/笔(哟)* pR-yPo *页)/笔(问:*(哟/笔(哟)* pR-yPo *页));
结束
结束

登录置评。

JK
JK 2018年2月8日
嗨Torsten,
谢谢你的时间花在这个问题上。整个问题是关于计算哟,滞留物流与边界值解算器x = 0。归来的成分渗透x = 0时,它的单位是摩尔/摩尔,而流动滞留物y是在摩尔/ s。我认为可能会让人困惑。历史上y在膜科学用于渗透成分,在普遍的气体组分热力学和因变量的数学。
ODE杭钢流膜,它是permability (Q)乘以部分压差/膜(压力retentae *摩尔分数在滞留物-压力渗透*渗透的摩尔分数)。滞留物年中的摩尔分数
年= nR /笔(nR) (nR y在我的程序——这是令人困惑的,我承认)
摩尔流渗透在x可以calulated平衡从x = 0到x
nP = nR-nRo
所以
yP = (nR-nRo) / (nR-nRo)和(或(y-yRo) / (y-yRo)和代码)
在x = 0 nr = nro,渗透流量是零和南交通等式变成了0/0。所以在x = 0我需要一个特殊的方程以避免NaN。在x = 0渗透流量为零,但渗透仍有一篇作文。它是
归来= n * /笔(n *)
n *是膜的流
n * = Q * (pR哟/ sum(哟)- pP归来)
这是解决常规PermOmega计算。
如果我能计算滞留物浓度直接在x = 0(哟),我就不会解决BV-problem。可能有一个代数解但当我添加的压力和温度资料,我需要解决ODE。
所以我怎么能改善bvp4c-routine的稳定?昨天我试着逐渐增加膜的数量和使用旧的解决方案(延续)作为一个新的起点。它没有工作。当我增加膜的数量从1.55到1.56系统转入混乱和光滑的轮廓将becomme混乱。我也改变了网格大小。更小的网格大小与更多指向x接近1会提高稳定性。
更多的想法如何使这项工作吗?它必须是可能的。
JK
JK 2018年2月9日
编辑:JK 2018年2月10日
今天晚上好,迫不及待地想试一试,会让你更新。
编辑:我得到的雅可比矩阵奇异解算器达到x = 1。在边界条件我不定义正确的边值哟。有趣的是我只能通过8个变量的边界函数。
我定义新的y (y_old;哟),所以新的y = y(1:4)和哟= y (8)。
这是代码,很简单:
清晰的
全球yalpha
yalpha = (0.3; 0.25; 0.01; 0.01);%进气摩尔流
n = 1;%没有膜
Q = [2.3369 e-08;1.17095 e-06;3.81087 e-08;2.36619 e-07] * n;%渗透率
yo0 = (0.95; 0.60; 0.92; 0.70)。* yalpha;%出口流量的估算值
%边值解算器
公元前bcfun = @(是的,yb)(是的,yb);
initguess = @ (x) init (x, yo0);
odefun = @ (x, y,哟)ydot (x, y, n);
solinit = bvpinit (linspace (0, 1, 11), initguess);
选择= bvpset (“统计数据”,“上”,“RelTol”1的军医);
索尔= bvp4c (odefun bcfun、solinit选项);
- - - - - - - - - - - - - %功能
函数IG = init (x, yo0)%网格的初始猜测
全球yalpha
IG = 0 (8, 1);
搞笑(1:4)= yo0 + x * (yalpha-yo0);
搞笑(8)= yo0;
结束
函数res = bc (ya,)%边值函数
全球yalpha
res = 0 (8, 1);
res(1:4) =丫(1:4)丫(8);%左BV (x = 0) y =哟
res (8) = yb (1:4) -yalpha;% y右侧BV (x = 1) y = yalpha
结束
函数就要= ydot (x, y, n)% '
就要= 0 (8,1);
Q = [2.3369 e-08;1.17095 e-06;3.81087 e-08;2.36619 e-07] * n;
公关= 1 e6;
页= 1 e5;
它们(1:4)=问:* (y(1:4) /笔(y(1:4)) *公关- (y (1:4) - y(8)) /笔(y (1:4) - y(8)) *页);
它们(8)= 0 (4,1);
结束
我想我们差不多了……
2的评论
显示1年长的评论
JK
JK 2018年2月12日
雅可比矩阵奇异。yb(8)不是边界函数中定义。服饰业/ dyb(8)线的雅可比矩阵边界值是空的。

登录置评。

JK
JK 2018年2月10日
恩:
问题的解决者归来x = 0。当n > > 1有时想出解决方案,min(归来)< 0或max(归来)> 金宝搏官方网站1 phyiscally毫无意义。所以我calcualte归来为n = 1或任何一个稳定的解决方案,然后使用它作为一个新的startingpoint未来计算的归来。我使用延续扩展解决方案所需的膜。
这是代码:
yalpha = (0.3; 0.25; 0.01; 0.01);摩尔%输入流
N = 5;%的膜
全球yPo0
Q = [2.3369 e-08;1.17095 e-06;3.81087 e-08;2.36619 e-07];渗透率
yPo0 = Q /笔(Q);最初的猜一猜
%值估计的膜(x = 0)
yo0 = (0.95; 0.60; 0.92; 0.70)。* yalpha;
%开始
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
抽搐;
%边值解算器
xx = linspace (0, 1, 10);
xx = xx。^ 0.5;%移动网格指向x = 1
图
BCJac ={(诊断接头(1(1、4));0(4,4)],[0(4,4);诊断接头((1、4))],[诊断接头(1 (1,4);0 (4,4)]};
bcfun = @(是的,yb,哟)[ya-yo;yb-yalpha];
initguess = @ (x)搞笑(x, yo0 yalpha);
索尔= bvpinit (xx, initguess yo0);
nn = (n - 1) / 10;
为n = 1:神经网络:n
%渗透率
Q = [2.3369 e-08;1.17095 e-06;3.81087 e-08;2.36619 e-07] * n;
FJac = @ (x, y,哟)Jacfun (x, y,哟,n);
odefun = @ (x, y,哟)ydot (x, y,哟,n);
选择= bvpset (“统计数据”,“上”,“RelTol”1 e - 3,“BCJacobian”,BCJac);%的FJacobian FJac,
索尔= bvp4c (odefun bcfun、溶胶、期权);
yy =德瓦尔(溶胶,xx);
情节(xx, yy (2:));
轴([0 1 0,0.3]);
标题(num2str (n));
暂停(0.2);
索尔= bvpinit (xx, @ (x)德瓦尔(溶胶,x), sol.parameters);%确保网不是膨胀与计算网格点
结束
disp (sol);
toc;
x = sol.x;
y = sol.y;
情节(x, y);
标题([num2str (n),“膜”]);
包含(“无量纲长度”);
ylabel (“摩尔流滞留物”);
函数就要= ydot (x, y,哟,n)
Q = [2.3369 e-08;1.17095 e-06;3.81087 e-08;2.36619 e-07] * n;
公关= 1 e6;
页= 1 e5;
如果x = = 0
全球yPo0
选择= optimset (“显示”,“没有”);
如果分钟(yPo0) < 0 | |马克斯(yPo0) 1
yPo0 = (0.2; 0.75; 0.03; 0.02);
哔哔的声音;%如果你听到:不好
结束
归来= fsolve(@(归来)PermOmega(归来),yPo0,选项);
yPo0 =归来;%的突破!
就要=问:* (y /笔(y) * pR-yPo *页);
其他的
就要=问:* (y /笔(y) *公关- (y-yo) / sum (y-yo) *页);
结束
函数f = PermOmega(归来)
%渗透成分x = 0
f = yPo-Q。*(哟/笔(哟)* pR-yPo *页)/笔(问:*(哟/笔(哟)* pR-yPo *页));
结束
结束
函数[FJ, FJo] = Jacfun (x, y,哟,n)
FJ =洁西(x, y,哟,n);
FJo = Jacyo (x, y,哟,n);%对雅克比外部函数
结束
函数f =搞笑(x, yo0 yalpha)%初始猜测只用于第一个计算
f = yo0 + x * (yalpha-yo0);
f (2) = yo0 (2) + x ^ 0.5 * (yalpha (2) -yo0 (2));%需要y的凸性质(2)考虑在内
结束
它是有趣的,他们告诉你的第一件事是,当数值计算解决方案是起始值。金宝搏官方网站和你认为你理解但打夜盯着电脑屏幕后你终于意识到:这一切都是开始值!
由于Torsten,你有很大的帮助!
纳迪姆阿巴斯
纳迪姆阿巴斯 2019年9月3日
你好亲爱的
如何找到BVP4C方法的解决方案的稳定性?金宝搏官方网站
1评论
的孩子叫K G
的孩子叫K G 2023年5月31日
根据文学作品通过不同的初始猜测

登录置评。

类别

找到更多的在方程求解在帮助中心和文件交换

标签

下载188bet金宝搏

社区寻宝

找到宝藏在MATLAB中央,发现社区如何帮助你!

开始狩猎!