用Matlab求解一个六阶非线性微分方程

24次浏览(过去30天)

显示旧的注释

Wissem-Eddine KHATLA 2022年4月6日

0
链接

直接链接到这个问题

https://au.mathworks.com/matlabcentral/answers/1689880-solving-a-6th-order-non-linear-differential-equation-in-matlab

评论道: Wissem-Eddine KHATLA2022年4月20日

大家好,

我正试图解一个高阶非线性微分方程

：

作为边界条件，我有以下几个条件:

我试着用一种射击法算法来解这个方程

，我把它转换成一阶微分方程的形式如下:

根据边界条件，有:

我们要猜剩下的几个:

我在我的帖子中附上了我写的算法，它受到传记研究的启发:主要代码称为“shoot4nl”。m”，我试着猜测一些初始值，以便运行算法。我还选择了一个值

非常“高”，以匹配条件……

通过这样做，我得到以下错误:

                          警告:矩阵为奇异的工作精度。
                         
                          在newtonRaphson2中(第23行)
                         
                          在射击4(第14行)
                         
                          错误使用newtonRaphson2
                         
                          太多次迭代
                         
                          错误在shoot4n1(第14行)中
                         
                          u = newtonRaphson2(@residual,u);

我认为这意味着在这个过程中会出现一个除以0的情况，但我仍然不知道如何解决这个问题……

有没有人能帮我一下，或者把他的数学专业知识带来给我一些提示?

提前谢谢你，

致以最亲切的问候。

2的评论
显示隐藏 1旧评论

Wissem-Eddine KHATLA 2022年4月7日

@Bruno陈德良谢谢你的回复。 @Torsten 通过使用各种内置函数为我的第一次尝试提供了更正，但似乎我得到了一个错误，说明“fsolve”函数不正确……

登录评论。

登录回答这个问题。

答案(2)

Torsten 2022年4月6日

1
链接

直接链接到这个答案

https://au.mathworks.com/matlabcentral/answers/1689880-solving-a-6th-order-non-linear-differential-equation-in-matlab#answer_936070

编辑:Torsten 2022年4月6日

我认为，如果有经过良好测试的MATLAB替代品，这个论坛中没有人会潜入自制软件。

所以这几行代码可以让你开始使用专业软件:

                              b0 = [1 1 1 1];
                             
                              BC = fsolve(@fun,bc0);
                             
                              F0 = [1 -5;0;bc(1);bc(2);bc(3);bc(4)]
                             
                              [-10,10];
                             
                              [eta,F] = 0 (@fun_ode,etaspan, 0);
                             
                              情节(η,F (: 1))
                             
                              函数Res =乐趣(bc)
                             
                              F0 = [1e-5;0;bc(1);bc(2);bc(3);
                             
                              Etaspan = [- 10,0,10];
                             
                              [eta,F] = 0 (@fun_ode,etaspan, 0);
                             
                              res(1) = F(2,2);
                             
                              res(2) = F(2,4);
                             
                              res(3) = F(3,1);
                             
                              res(4) = F(3,2);
                             
                              结束
                             
                              函数dfdelta = fun_ode(eta,F)
                             
                              dFdeta = [F (2), (3); F (4), F (5), F (6); (1/9 * (5 * F(1) 4 *η* F (2)) 3 * F (1) ^ 2 * (2) * F (6)) / F (1) ^ 3];
                             
                              结束

30的评论
显示隐藏 29旧评论

Torsten 2022年4月6日

编辑:Torsten 2022年4月6日

下面的代码可以工作，例如:，和gives the included solution for F:

                                   Bc0 = [2.2504e+06 -2.5004e+05 -2.2123e+06 5.6616e+05]
                                  
                                   Etaspan = [-30 30];
                                  
                                   BC = fsolve(@(BC)fun(BC,etaspan)，bc0);
                                  
                                   bc(1);bc(2);bc(3);
                                  
                                   选项= odeset()“InitialStep”1 e-9);%, ' MaxOrder ', 1);
                                  
                                   [eta,F] = ode45(@fun_ode,etaspan,bc_total);%,选项);
                                  
                                   情节(η,F (: 1))
                                  
                                   函数Res = fun(bc,etaspan)
                                  
                                   bc(1);bc(2);bc(3);
                                  
                                   etasp_ode = [etaspan(1)，0,etaspan(2)];
                                  
                                   选项= odeset()“InitialStep”1 e-9);%, ' MaxOrder ', 1);
                                  
                                   [eta,F] = ode45(@fun_ode,etaspan_ode,bc_total);%,选项);
                                  
                                   res(1) = F(2,2);
                                  
                                   res(2) = F(2,4);
                                  
                                   res(3) = F(3,1)-1e-7;
                                  
                                   res(4) = F(3,2);
                                  
                                   结束
                                  
                                   函数dfdelta = fun_ode(eta,F)
                                  
                                   dFdeta = [F (2), (3); F (4), F (5), F (6); (1/9 * (5 * F(1) 4 *η* F (2)) 3 * F (1) ^ 2 * (2) * F (6)) / F (1) ^ 3];
                                  
                                   结束

Torsten 2022年4月7日

使用下面的代码，将其命名为“main”。M”，在八度编辑器中加载并运行它。

                                   函数主要
                                  
                                   Bc0 = [2.2504e+06 -2.5004e+05 -2.2123e+06 5.6616e+05]
                                  
                                   Etaspan = [-30 30];
                                  
                                   BC = fsolve(@(BC)fun(BC,etaspan)，bc0);
                                  
                                   bc(1);bc(2);bc(3);
                                  
                                   选项= odeset()“InitialStep”1 e-9);%, ' MaxOrder ', 1);
                                  
                                   [eta,F] = ode45(@fun_ode,etaspan,bc_total);%,选项);
                                  
                                   情节(η,F (: 1))
                                  
                                   结束
                                  
                                   函数Res = fun(bc,etaspan)
                                  
                                   bc(1);bc(2);bc(3);
                                  
                                   etasp_ode = [etaspan(1)，0,etaspan(2)];
                                  
                                   选项= odeset()“InitialStep”1 e-9);%, ' MaxOrder ', 1);
                                  
                                   [eta,F] = ode45(@fun_ode,etaspan_ode,bc_total);%,选项);
                                  
                                   res(1) = F(2,2);
                                  
                                   res(2) = F(2,4);
                                  
                                   res(3) = F(3,1)-1e-7;
                                  
                                   res(4) = F(3,2);
                                  
                                   结束
                                  
                                   函数dfdelta = fun_ode(eta,F)
                                  
                                   dFdeta = [F (2), (3); F (4), F (5), F (6); (1/9 * (5 * F(1) 4 *η* F (2)) 3 * F (1) ^ 2 * (2) * F (6)) / F (1) ^ 3];
                                  
                                   结束

Torsten 2022年4月8日

编辑:Torsten 2022年4月8日

为什么选择在函数中计算4个残差?

正如我已经试图向您解释的那样，我从eta = -eta0 (F = F' = 0)的两个条件开始ode积分器，并对F''，F''，F''''和F'''''的剩余四个初始条件进行估计。你定义了另外四个条件(两个在0处，两个在= eta0处)应该成立。所以我们有四个值需要消去，即F' = F'' =0在=0时F = F' =0在=eta0时。这是需要求解的四个方程。在res向量中，错误被返回到fsolve，并提示为F " (-eta0)，…，F'''''(-eta0)提供更好的猜测，以消除这些错误。

我明白这句话的目的:

BC = fsolve(@(BC)fun(BC,etaspan)，bc0);

然而，我不熟悉这种语法(匿名函数):为什么你选择这样写?

“fun”中需要“etaspan”来定义ode积分器的积分间隔。这条线

BC = fsolve(@(BC)fun(BC,etaspan)，bc0);

具有将“etaspan”作为可在函数中使用的附加参数传递给“fun”的效果。如果你不喜欢这样，你将不得不在调用函数和“fun”中同时更改“etaspan”，如果你更改了eta0。

我非常希望“bvp4c”能比我提供的ad-hoc拍摄方法更稳定地解决你的问题。你已经在MATLAB中测试过我给你的代码了吗?

Torsten 2022年4月8日

编辑:Torsten 2022年4月8日

由于我没有MATLAB可用，我无法测试。

但是这个设置(eta0 = 30, F(eta0)=1e-7)对于求解者来说已经很困难了。从不那么激进的值开始(eta0 = 10, F(eta0) = 1e-3或其他)。如果成功，使用该解作为初始值，使eta0变大，F(eta0)变小。

但同样在OCTAVE下，计算非常不稳定，eta0和F(eta0)的微小变化会导致eta=0时高度的巨大变化。虽然您得到了奇异雅可比矩阵的错误消息，但您应该使用我给您的代码尝试bvp4c或bvp5c。一般来说，用有限差分法求解边值问题要比射击稳定得多。

在有限区间内求解(如[-1:1]，通过坐标变换eta~ = tanh(eta))可以稳定计算。

不知怎么的，我觉得问题的表述中缺少了一个条件。如果从x=0开始积分，有三个条件(因为我认为解与=0对称)即F'(0)=F'' (0)=F'''''(0)=0。另外，在无穷远处有两个条件:F(无穷)= F'(无穷)= 0。但第六个必要条件是什么?

Torsten 2022年4月9日

编辑:Torsten 2022年4月9日

你认为这是一个很好的数值公式吗?

你知道这是一个完全不同的问题吗?旧公式的结果在eta = 0时给出了约为1e6的值。现在你规定F = 1在= 0时。你认为这两个公式除了基本的微分方程之外还有什么共同之处吗?

这应该是重新表述问题的正确代码，但它不收敛:

                                   Bc0 = [-2.5004e+05 -2.2123e+06 5.6616e+05]
                                  
                                   Etaspan = [0 10];
                                  
                                   BC = fsolve(@(BC)fun(BC,etaspan)，bc0);%在eta = eta0时调整F''，F''''，F'''''的初始条件
                                  
                                   b_total = [1;0;bc(1);0;bc(2);
                                  
                                   [eta,F] = ode45(@fun_ode,etaspan,bc_total);
                                  
                                   情节(η,F (: 1))
                                  
                                   %为拍摄方法创建残差
                                  
                                   函数Res = fun(bc,etaspan)
                                  
                                   b_total = [1;0;bc(1);0;bc(2);
                                  
                                   etasp_ode = [etaspan(1)，etaspan(2)];
                                  
                                   [eta,F] = ode45(@fun_ode,etaspan_ode,bc_total);
                                  
                                   res(1) = F(end,1)-1e-2;%在eta = eta0时F与0的偏差
                                  
                                   res(2) = F(end,2);%在eta = eta0时F'与0的偏差
                                  
                                   res(3) = F(end,3);在eta = eta0时F''与0的偏差%
                                  
                                   结束
                                  
                                   % 1- 1阶ODE系统
                                  
                                   函数dfdelta = fun_ode(eta,F)
                                  
                                   dFdeta = [F (2), (3); F (4), F (5), F (6); (1/9 * (5 * F(1) 4 *η* F (2)) 3 * F (1) ^ 2 * (2) * F (6)) / F (1) ^ 3];
                                  
                                   结束

Wissem-Eddine KHATLA 2022年4月9日

@Torsten ，

我在之前的帖子中犯了一个错误:我正在听从你的建议，试图解决问题

只间隔。在这种情况下，我引入了一个新的边界条件，如:

我认为这是解决问题的唯一方法，因为正如你所注意到的:用OCTAVE获得的解决方案真的不稳定。在= 0处的解值非常重要，所以我需要一个可靠的方法来确定它。

以下是对脚本的修改:

                                   %主脚本的拍摄方法
                                  
                                   b0 = [2.2504e+06 2.212e3 -6 5.6616e-05];
                                  
                                   Etaspan = [0 10];
                                  
                                   BC = fsolve(@(BC)fun(BC,etaspan)，bc0);%在eta = eta0时调整F''，F''''，F'''''的初始条件
                                  
                                   b_total = [bc(1);0;bc(2);0;
                                  
                                   [eta,F] = ode45(@fun_ode,etaspan,bc_total);
                                  
                                   情节(η,F (: 1))
                                  
                                   %为拍摄方法创建残差
                                  
                                   函数Res = fun(bc,etaspan)
                                  
                                   b_total = [bc(1);0;bc(2);0;
                                  
                                   etasp_ode = [etaspan(1)，etaspan(2)];
                                  
                                   [eta,F] = ode45(@fun_ode,etaspan_ode,bc_total);
                                  
                                   res(1) = F(end,1)-1e-2;%在eta = eta0时F与0的偏差
                                  
                                   res(2) = F(end,2);%在eta = eta0时F'与0的偏差
                                  
                                   res(3) = F(end,3);在eta = eta0时F''与0的偏差%
                                  
                                   结束
                                  
                                   % 1- 1阶ODE系统
                                  
                                   函数dfdelta = fun_ode(eta,F)
                                  
                                   dFdeta = [F (2), (3); F (4), F (5), F (6); (1/9 * (5 * F(1) 4 *η* F (2)) 3 * F (1) ^ 2 * (2) * F (6)) / F (1) ^ 3];
                                  
                                   结束

总是有一个错误……

我还要试试我以前的配方

)，但

这是一个很难解决的问题……

再次感谢你。

Torsten 2022年4月9日

编辑:Torsten 2022年4月9日

在= 0处的解值非常重要，所以我需要一个可靠的方法来确定它。

我不想打击您的信心，但是对旧代码的实验表明，eta=0处的值与所选择的eta0有很大的关系。

这在OCTAVE下工作:

                                   b =[3.8309 0.4494 -0.1066];
                                  
                                   Etaspan = [0 12];
                                  
                                   BC = fsolve(@(BC)fun(BC,etaspan)，bc0);%在eta = eta0时调整F''，F''''，F'''''的初始条件
                                  
                                   b_total = [bc(1);0;bc(2);0;
                                  
                                   选项= odeset()“InitialStep”1 e-9);
                                  
                                   [eta,F] = ode45(@fun_ode,etaspan,bc_total,options);
                                  
                                   公元前
                                  
                                   情节(η,F (: 1))
                                  
                                   结束
                                  
                                   %为拍摄方法创建残差
                                  
                                   函数Res = fun(bc,etaspan)
                                  
                                   b_total = [bc(1);0;bc(2);0;
                                  
                                   etasp_ode = [etaspan(1)，etaspan(2)];
                                  
                                   选项= odeset()“InitialStep”1 e-9);
                                  
                                   [eta,F] = ode45(@fun_ode,etaspan_ode,bc_total,options);
                                  
                                   res(1) = F(end,1)-1e-5;%在eta = eta0时F与0的偏差
                                  
                                   res(2) = F(end,2);%在eta = eta0时F'与0的偏差
                                  
                                   res(3) = F(end,3);在eta = eta0时F''与0的偏差%
                                  
                                   结束
                                  
                                   % 1- 1阶ODE系统
                                  
                                   函数dfdelta = fun_ode(eta,F)
                                  
                                   % F
                                  
                                   dFdeta = [F (2), (3); F (4), F (5), F (6); (1/9 * (5 * F(1) 4 *η* F (2)) 3 * F (1) ^ 2 * (2) * F (6)) / F (1) ^ 3];
                                  
                                   结束

Wissem-Eddine KHATLA 2022年4月10日

编辑:Wissem-Eddine KHATLA 2022年4月10日

@Torsten 我尝试了bvp4c求解器使用拍摄方法的结果:它似乎不能很好地使用以下代码:

                                   %求解方法采用bvp4c求解器
                                  
                                   %边界条件:F (0) = F (0) = F”的(0)= 0 (eta0) = F (eta0) = F”(eta0) = 0
                                  
                                   Eta0 = 20;
                                  
                                   Eta = linspace(0,eta0,100);
                                  
                                   Yinit = [3.3517;0;0.5079;0;0.1194;0];OCTAVE ode45集成的%解决方案
                                  
                                   Solinit = bvpinit(eta, init);
                                  
                                   Sol = bvp4c(@fun_ode， @bcfun, solinit);
                                  
                                   Eta = sol.x;
                                  
                                   F = sol.y(1，:);
                                  
                                   数字
                                  
                                   情节(埃塔(F);
                                  
                                   函数dF = fun_ode(eta,F)
                                  
                                   dF = [F (2);
                                  
                                   F (3);
                                  
                                   F (4);
                                  
                                   F (5);
                                  
                                   F (6);
                                  
                                   (1/9 *(5 *(1) 4 *η* F (2)) 3 * F (1) ^ 2 * (2) * F (6)) / F (1) ^ 3);
                                  
                                   结束
                                  
                                   函数Res = bcfun(ya,yb)
                                  
                                   res =[丫(2);丫(4);丫(6);yb (1);yb (2);yb (3)];
                                  
                                   结束

另一位用户建议我这个问题并不棘手:我很困惑，因为文献中有一个类似的方程，所以我认为我仍然有可能解决它。

谢谢你！

登录评论。

布鲁诺陈德良 2022年4月9日

0
链接

直接链接到这个答案

https://au.mathworks.com/matlabcentral/answers/1689880-solving-a-6th-order-non-linear-differential-equation-in-matlab#answer_938920

使用bvp4c

Eta0 = 10;

Eta = linspace(0,eta0,100);

[0;0;0;0];

Solinit = bvpinit(eta, init);

Sol = bvp4c(@fun_ode， @bcfun, solinit);

Eta = sol.x;

F = sol.y(1，:);

情节(埃塔(F);

函数dF = fun_ode(eta,F)

dF = [F (2);

F (3);

F (4);

F (5);

F (6);

(1/9 *(5 *(1) 4 *η* F (2)) 3 * F (1) ^ 2 * (2) * F (6)) / F (1) ^ 3);

结束

函数Res = bcfun(ya,yb)

res =[丫(2);% F (0)

丫(4);% F”(0)

丫(6);% F”“(0)

yb (1:3)];% F(eta0)， F'(eta0)， F'(eta0)，

结束

79条评论
显示隐藏 78旧评论

布鲁诺陈德良 2022年4月10日

“然而，你最后的绘图提供了良好的形状(良好的值?)。”

我不确定你说的“外形好”和“性价比高”是什么意思。我不知道解析解，也不知道其他5个导数值是如何正确匹配的。对五阶导数施加条件简直是疯了。导数是一个非常不稳定的算子，取它5 / 6次，对它进行计算是没有希望的。

“选择‘yinit’有很强的依赖性……

因为非线性，我已经告诉过你们了。

当更多的算法告诉雅可比矩阵是奇异的，它只是表明问题是病态的，最终状态(或至少它的一些子空间)对初始状态的依赖在数字上是零。这可能是由于您的问题是严重的规模(对于大eta0 >> 1，正如我上面所说的)或在您试图解决的问题的本质。

另外，你能告诉我为什么你要这样写残差吗?

Res =[ya([2 4 6]);yb (1:3)];

从bc你告诉我们，如果你错过了我的评论在这里的第一个代码

                                   res =[丫(2);% F (0) = 0
                                  
                                   丫(4);% F”(0)= 0
                                  
                                   丫(6);% F”“(0)= 0
                                  
                                   yb (1:3)];% F(eta0)=0, F'(eta0)=0, F'(eta0)=0，

Wissem-Eddine KHATLA 2022年4月11日

编辑:Wissem-Eddine KHATLA 2022年4月11日

@Bruno陈德良这是从所谓的“润滑理论”中描述高度场的方程推导出来的:高阶是因为我们考虑了一个4阶的压力分量。

物理问题涉及间隔上的6 BC

它们是:

@Torsten 前面的代码给出了一个类似于预期的形状，但我不知道如何为我的问题添加新的约束，以减少最终结果与初始猜测或值之间的依赖关系 η ．

                                   %求解方法采用bvp4c求解器
                                  
                                   %边界条件:F (0) = F (0) = F”的(0)= 0 (eta0) = F (eta0) = F”(eta0) = 0
                                  
                                   Eta0 = 40;
                                  
                                   [1;0;0.5079;0;0.1194;0];bvpinit解算器的初始条件。
                                  
                                   Nbiter = 10;%迭代次数。
                                  
                                   我们选择以对数尺度迭代eta的值
                                  
                                   我们使用bvp4c求解器对相应的ODE进行积分
                                  
                                   eta0 = logspace(0,log10(eta0)，nbiter);
                                  
                                   为k = 1: nbiter
                                  
                                   Eta = linspace(0,eta0tab(k)，100);
                                  
                                   Solinit = bvpinit(eta, init);
                                  
                                   Sol = bvp4c(@fun_ode， @bcfun, solinit);
                                  
                                   Eta = sol.x;
                                  
                                   F = sol.y(1，:);
                                  
                                   yit = sol.y(:，1);
                                  
                                   结束
                                  
                                   我们绘制结果
                                  
                                   情节(η,- f);
                                  
                                   我们构造一个1阶ODE的向量
                                  
                                   函数dF = fun_ode(eta,F)
                                  
                                   dF = [F (2:6);
                                  
                                   (1/9 *(5 *(1) 4 *η* F (2)) 3 * F (1) ^ 2 * (2) * F (6)) / F (1) ^ 3);
                                  
                                   结束
                                  
                                   我们用y(b)-a = 0建立该方法的残差。
                                  
                                   函数Res = bcfun(ya,yb)
                                  
                                   Res =[ya([2 4 6]);yb (1:3)];
                                  
                                   结束

谢谢你的帮助。

Wissem-Eddine KHATLA 2022年4月15日

@Bruno陈德良， @Torsten 谢谢你的贡献。我有了一个新的想法并且对原来的方程做了一些修改。我介绍了:

实验常数一个
约束:施加在体积上的约束，产生一个新的方程:

这两个一个和问是已知的。因此，我得到这两个方程作为ODE的系统来实现:

我只是在最后一行中将第二个方程添加到ODE的系统中

                                   dF = [F (2:6);
                                  
                                   (1 / (9 *)) * (5 * F(1) 4 *η* F (2)) 3 * F (1) ^ 2 * (2) * F (6)) / F (1) ^ 3;F (1)];

边界条件仍然相同:

我只是在积分上加上一个条件

我想这就得到了一个适定的问题。

如果我测试修改后的代码:A = Q = 1:我得到一个绘图。然而，如果我选择应用A和Q的真实值:代码再次崩溃……

                                   在[0;eta0]区间内求解弹性区域内弹性液滴的扩散问题
                                  
                                   %求解方法采用bvp4c求解器
                                  
                                   %边界条件:F (0) = F (0) = F”的(0)= 0 (eta0) = F (eta0) = F”(eta0) = 0
                                  
                                   %对卷的约束
                                  
                                   Eta0 = 50;
                                  
                                   [1];0;0.5079;0;0.1194;0;1);bvpinit解算器的初始条件。
                                  
                                   %yinit = [0.25;0;1;0;0.5;0;1];
                                  
                                   Nbiter = 10;%迭代次数。
                                  
                                   我们选择以对数尺度迭代eta的值
                                  
                                   我们使用bvp4c求解器对相应的ODE进行积分
                                  
                                   eta0 = logspace(0,log10(eta0)，nbiter);
                                  
                                   为k = 1: nbiter
                                  
                                   Eta = linspace(0,eta0tab(k)，100);
                                  
                                   Solinit = bvpinit(eta, init);
                                  
                                   Sol = bvp5c(@fun_ode， @bcfun, solinit);
                                  
                                   Eta = sol.x;
                                  
                                   F = sol.y(1，:);
                                  
                                   yit = sol.y(:，1);
                                  
                                   结束
                                  
                                   %我们绘制结果[F;[F]在一个新的数字
                                  
                                   图()
                                  
                                   情节(η,- f,“- - -”，“线宽”2);
                                  
                                   网格在
                                  
                                   持有在
                                  
                                   情节(η,sol.y (2:)“线宽”, 2)
                                  
                                   传奇(“F”，“F”，“位置”，“东北”）
                                  
                                   包含(“\埃塔”，“字形大小”13);
                                  
                                   ylabel (“(\η),F (\ eta)”，“字形大小”, 13)
                                  
                                   我们首先构造一个1阶ODE的向量
                                  
                                   最后一个是约束项
                                  
                                   函数dF = fun_ode(eta,F)
                                  
                                   A = 1;% 3.7926 e;实验常数
                                  
                                   dF = [F (2:6);
                                  
                                   (1 / (9 *)) * (5 * F(1) 4 *η* F (2)) 3 * F (1) ^ 2 * (2) * F (6)) / F (1) ^ 3;F (1)];
                                  
                                   结束
                                  
                                   我们用y(b)-a = 0建立该方法的残差。
                                  
                                   函数Res = bcfun(ya,yb)
                                  
                                   Q = 1;% 0.4 e-9;%实验流速(m^3/s)
                                  
                                   Res =[ya([2 4 6]);yb (1:3);yb (7) - q);
                                  
                                   结束

你们对这个问题有什么想法吗?我的修改写对了吗?还是我的配方又出了问题?

再次感谢您的帮助 @Torsten 和 @Bruno陈德良

布鲁诺陈德良 2022年4月15日

编辑:布鲁诺陈德良 2022年4月15日

这段代码不会因为你的实验数据而崩溃，但是有一个警告。实际上你的系统是6阶的(假装是7阶的)你想施加6bc + 1积分约束，它不能满足所有的条件。

                                   清晰的
                                  
                                   关闭所有
                                  
                                   在[0;eta0]区间内求解弹性区域内弹性液滴的扩散问题
                                  
                                   %求解方法采用bvp4c求解器
                                  
                                   %边界条件:F (0) = F (0) = F”的(0)= 0 (eta0) = F (eta0) = F”(eta0) = 0
                                  
                                   %对卷的约束
                                  
                                   Eta0 = 50;
                                  
                                   A = 3.7926e-7;
                                  
                                   Q = 0.4e-9;
                                  
                                   %Finit = [1;0;0.5079;0;0.1194;0;1);bvpinit解算器的初始条件。
                                  
                                   有限扩散= [-1.0513 e + 04; 0, 158.4915, 0; -3.6654; 0; -1.0513 e + 04];
                                  
                                   Nbiter = 20;%迭代次数。
                                  
                                   我们选择以对数尺度迭代eta的值
                                  
                                   我们使用bvp4c求解器对相应的ODE进行积分
                                  
                                   Eta0tab = logspace(log10(1)，log10(eta0)， nbiter);
                                  
                                   为k = 1: nbiter
                                  
                                   流(“迭代% d / % d \ n”， k，苦);
                                  
                                   etascale = A^(1/6)
                                  
                                   Eta = linspace(0,eta0tab(k)，100);
                                  
                                   Ginit = Finit .* etascale。^ (0:5 0) ';
                                  
                                   solinit = bvpinit(xi,Ginit);
                                  
                                   试一试
                                  
                                   sol = bvp5c(@(xi, G) fun_ode(xi,G,A,Q,etascale)，．..
                                  
                                   @(是的,yb) bcfun(是的,yb, Q etascale),．..
                                  
                                   solinit);
                                  
                                   抓
                                  
                                   %如果k==nbiter
                                  
                                   % fprintf('最后一次迭代失败\n');
                                  
                                   %返回
                                  
                                   %结束
                                  
                                   继续
                                  
                                   结束
                                  
                                   Finit = sol.y(:，1) ./ etascale. /^ (0:5 0) ';
                                  
                                   结束
                                  
                                   Eta = etascale*xi;
                                  
                                   F = sol.y(1，:);
                                  
                                   Fp = sol.y(2，:) / etascale;
                                  
                                   %我们绘制结果[F;[F]在一个新的数字
                                  
                                   图()
                                  
                                   情节(η,- f,“- - -”，“线宽”2);
                                  
                                   网格在
                                  
                                   持有在
                                  
                                   情节(η,《外交政策》,“线宽”, 2)
                                  
                                   传奇(“F”，“F”，“位置”，“东北”）
                                  
                                   包含(“\埃塔”，“字形大小”13);
                                  
                                   ylabel (“(\η),F (\ eta)”，“字形大小”, 13)
                                  
                                   我们首先构造一个1阶ODE的向量
                                  
                                   最后一个是约束项
                                  
                                   函数dG = fun_ode(xi,G,A,Q,etascale)
                                  
                                   Eta = xi*etascale;
                                  
                                   F = G ./ (etascale)。^ (0:5 0) ';
                                  
                                   dF = [F (2:6);
                                  
                                   (1 / (9 *)) * (5 * F(1) 4 *η* F (2)) 3 * F (1) ^ 2 * (2) * F (6)) / F (1) ^ 3;
                                  
                                   F (1) /);
                                  
                                   dG = dF .* (etascale)。1 ^[1:6]”;
                                  
                                   结束
                                  
                                   我们用y(b)-a = 0建立该方法的残差。
                                  
                                   函数res = bcfun(ya,yb,A,Q,etascale)
                                  
                                   Ya = Ya ./ (etascale)。^ (0:5 0) ';
                                  
                                   Yb = Yb ./ (etascale)。^ (0:5 0) ';
                                  
                                   Res =[ya([2 4 6]);yb (1:3);yb (7) - q / A];
                                  
                                   结束

布鲁诺陈德良 2022年4月16日

编辑:布鲁诺陈德良 2022年4月16日

@Wissem-Eddine KHATLA 这是不正确的匿名函数语法。

@guess(η,eta0 Q)

看起来你想写(没有@)

bvpinit(η,猜(η,eta0 Q))

这是向量的初始化语法。

IMO Torsen想用函数 bcpinit 因为我不会在跨度上复制值医生

函数-对于给定的网格点，guess函数必须返回一个向量，其元素是对解的相应组件的猜测。函数必须是这样的形式Y = guess(x)x是一个网格点和吗y是一个向量，其长度与解中分量的个数相同。例如，如果yinit是一个函数，那么在每个网格点bvpinit调用Y (:，j) = guess(x(j))对于多点边值问题，猜测函数必须为Y = guess(x, k)y对解的初步猜测是x在地区k．函数必须接受输入参数k，为编写guess函数提供了灵活性。但是，该功能不是必须使用的k．

我已经尝试了所有Torsen的想法，而不是不贴在这里，因为他们没有带来任何东西，不是说它降低了稳定性。

在玩了一段时间之后，我认为强迫积分是个错误的好主意。

Wissem-Eddine KHATLA 2022年4月20日

编辑:Wissem-Eddine KHATLA 2022年4月20日

@Torsten 是的，这是我的案子，但看起来他们融合在一起了

在本文中，我已经详细说明的条件。不幸的是，我的代码不工作，因为我应该找到

．..

代码如下:

                                   Eta0 = 20;
                                  
                                   Etamesh = [linspace(-eta0,0,100)，linspace(0,eta0,100)];
                                  
                                   Finit = [1;0.5;1;0;0);
                                  
                                   solinit = bvpinit(etamesh,Finit)
                                  
                                   Sol = bvp5c(@fun_ode， @bc, solinit);
                                  
                                   图()
                                  
                                   情节(sol.x sol.y (1:)“线宽”, 1.5)%绘图F
                                  
                                   ylim (2 [2]);
                                  
                                   持有在
                                  
                                   情节(sol.x sol.y (3:)“线宽”, 1.5)%绘图F”
                                  
                                   情节(sol.x sol.y (5:)“线宽”, 1.5)%绘图F''''
                                  
                                   网格在
                                  
                                   传奇(“F”，“F”“，“F”“”，“位置”，“西北”）
                                  
                                   函数dfdelta = fun_ode(eta,F,region)
                                  
                                   dFdeta = [F (2), (3); F (4), F(5),(行进(1))/ F (1) ^ 3];
                                  
                                   结束
                                  
                                   %边界条件:
                                  
                                   % F " '(-eta0) = F''''(-eta0) = 0
                                  
                                   % f (0) = 1
                                  
                                   % [1] = [2] = 0 on eta = eta0
                                  
                                   函数res = bc(FL,FR)
                                  
                                   = cos(4* /5) + i*sin(4* /5);
                                  
                                   Beta = conj(alpha);
                                  
                                   Gamma = abs(alpha);
                                  
                                   res = [FL(4,1);FL(5、1);．..
                                  
                                   Fr (1,1) - 1;．..
                                  
                                   FR (1, 1) fl(1、2);FR (2, 1) fl (2, 2);FR (3,1) fl (3 2);FR (4,1) fl (4,2);FR(5、1)fl (5, 2);．..
                                  
                                   (1 -(β+α))* FR(3 2) +(-(α+β)+γ^ 2)* FR (2, 2);(-(β+α)+γ^ 2)* FR (3 2) + (FR(2, 2) 1) *伽马^ 2);
                                  
                                   结束

我认为我的代码一切正常:也许你会发现一个错误?

谢谢你！

Wissem-Eddine KHATLA 2022年4月20日

@Torsten F''(-eta0)不等于1.35。例如，下面是我的情节:

与文章中的一个相比:

登录评论。

登录回答这个问题。

类别

MATLAB 数学数值积分与微分方程边值问题

社区寻宝

在MATLAB Central中寻找宝藏，并了解社区如何帮助您!

开始狩猎!

用Matlab求解一个六阶非线性微分方程

2的评论
显示隐藏 1旧评论

答案(2)

30的评论
显示隐藏 29旧评论

79条评论
显示隐藏 78旧评论

另请参阅

类别

标签

社区寻宝

用Matlab求解一个六阶非线性微分方程

2的评论 显示隐藏 1旧评论

答案(2)

30的评论 显示隐藏 29旧评论

79条评论 显示隐藏 78旧评论

2的评论
显示隐藏 1旧评论

30的评论
显示隐藏 29旧评论

79条评论
显示隐藏 78旧评论