求解一个参数未知的非线性ODE

7次浏览(过去30天)

显示旧的注释

khaoula Oueslati 2022年4月14日

0
链接

直接链接到这个问题

https://ch.mathworks.com/matlabcentral/answers/1696395-solving-a-non-linear-ode-with-unknown-parameter

评论道: khaoula Oueslati2022年4月19日

答:接受 Torsten

你好!我正在解决一个动力学参数未知A的ODE方程，我一直在使用python和深度学习来求解方程，并确定A的值，然而损失函数总是在10**4的数量级，参数A是错误的，我尝试了不同的超参数，但它不起作用。这是ODE方程: 顺铂/ dt = k1 * ((DP) ^ 2) 和k1 = k = 一个 e ^ ( 1 / R (- E / ( T + 273 ）） )， A是10**8的数量级，我有DP(t)数据。

我被困住了，我想知道用matlab解决这个问题的最佳方法是什么?或者还有其他类似的例子吗?

任何帮助都非常感谢!

0评论
显示隐藏-1旧的注释

登录评论。

登录回答这个问题。

接受的答案

Torsten 2022年4月19日

1
链接

直接链接到这个答案

https://ch.mathworks.com/matlabcentral/answers/1696395-solving-a-non-linear-ode-with-unknown-parameter#answer_945725

编辑:Torsten 2022年4月19日

                              %的时间点
                             
                              Ts =[1 2 3 4 5 6 7 8];
                             
                              Dp =[1000 700.32 580.42 408.20 317.38 281.18 198.15 100.12];
                             
                              P0 = 1e1;
                             
                              p = fminunc(@(p)fun(p,ts,DP)，p0)
                             
                              E = 111e3;
                             
                              R = 8.314;
                             
                              T = 371;
                             
                              A = p*exp(E/(R*T))
                             
                              情节(ts, DP)
                             
                              持有在
                             
                              情节(ts, 1. / (1 / DP(1) +一个* exp (e / (R * T)) * (ts-ts (1))));
                             
                              函数obj = fun(p,ts,DP)
                             
                              DP_model = 1/ (1/DP(1)+ p*(ts-ts(1)));
                             
                              obj = sum((DP-DP_model).^2)
                             
                              结束

6个评论
显示隐藏旧的评论

khaoula Oueslati 2022年4月19日

谢谢你！ @Torsten 感谢您的回答和宝贵的帮助!我真的很感激。

我尝试了你的代码，它给了我和你一样的输出，我尝试了另一个数据集，我知道A值，这是代码:

                                   回归数据百分比
                                  
                                   Dp = [1000 995.6179 991.2741 986.968 982.6991 978.467 974.2712 970.1113 965.9867 961.897 957.8419 953.8207 949.8332 945.8789 941.9574 938.0683 934.2111 930.3855 926.5912 922.8276 919.0945，.．.
                                  
                                   915.3915 911.7182 908.0743 904.4594 900.8732 897.3152 893.7853 890.283 886.8081 883.3602 879.939 876.5442 873.1755 869.8326 866.5152，.．.
                                  
                                   863.223 859.9557 856.7131 853.4948 850.3006 847.1302 843.9834 840.8599 837.7594 834.6817 831.6265，.．.
                                  
                                   828.5936 825.5827 822.5937 819.6262 816.68 813.755 810.8508 807.9673 805.1042 802.2614 799.4385 796.6354 793.852 791.0879，.．.
                                  
                                   788.343 785.6171 782.9099 780.2214 777.5512 774.8993 772.2654 769.6493 767.0509 764.47 761.9064 759.36 756.8305 754.3178，.．.
                                  
                                   751.8217 749.3421 746.8788 744.4316 742.0005 739.5851 737.1854 734.8013 732.4325 730.0789 727.7405 725.4169 723.1082，.．.
                                  
                                   720.814 718.5344 716.2692 714.0182 711.7813 709.5584 707.3494 705.154 702.9722 700.8039 698.6489 696.5072 694.3785 692.2628，.．.
                                  
                                   690.1599 688.0698 685.9923 683.9273 681.8748 679.8345 677.8063 675.7902 673.7861 671.7939 669.8134 667.8445 665.8872，.．.
                                  
                                   663.9413 662.0067 660.0834 658.1713 656.2701 654.38 652.5007 650.6321 648.7742 646.9269 645.0901 643.2637 641.4476];
                                  
                                   ts = linspace (1128128);
                                  
                                   E = 111e3;
                                  
                                   R = 8.314;
                                  
                                   T = 371;
                                  
                                   % ts=[1 2 3 4 5 6 7 8];
                                  
                                   % dp =[1000 700.32 580.42 408.20 317.38 281.18 198.15 100.12];
                                  
                                   情节(ts, DP)
                                  
                                   ylabel (“迪拜”）
                                  
                                   包含(“t”）
                                  
                                   持有在
                                  
                                   P0 = 1e1;
                                  
                                   p = fminunc(@(p)fun(p,ts,DP)，p0);
                                  
                                   disp ([“损失值为”num2str(乐趣(p, ts, DP))))
                                  
                                   disp ([“A值是”num2str (A)])
                                  
                                   情节(ts, 1. / (1 / DP (1) + p * (ts-ts (1))),“波”）;
                                  
                                   %我正在寻找的价值
                                  
                                   A = p*exp(E/(R*T));
                                  
                                   %目标函数
                                  
                                   函数obj = fun(p,ts,DP)
                                  
                                   DP_model = 1/ (1/DP(1)+ p*(ts-ts(1)));
                                  
                                   obj = sum((DP-DP_model).^2);
                                  
                                   结束

损失值为:损失值为0.35787， ground_truth的值为7.8e8

我不确定这种不匹配的来源。

顺便说一句，你是怎么找到DP_model表达式的?是某种近似吗?或者积分后我们得到解的表达式?

khaoula Oueslati 2022年4月19日

非常感谢 @Torsten ！并不是我不相信这个表达，我只是想知道如何找到它，这样我就可以向我的学术导师解释它。关于不匹配，可能是由于数据或E值的一些变化，我会进一步调查，再次感谢!

登录评论。

更多答案(3)

Torsten 2022年4月14日

1
链接

直接链接到这个答案

https://ch.mathworks.com/matlabcentral/answers/1696395-solving-a-non-linear-ode-with-unknown-parameter#answer_943080

D_p的ODE

(1/D_p0 + k1*(t-t0))

其中，D_p0 = D_p(t0)。

现在您可以应用“lsqcurvefit”来拟合未知参数A。

2的评论
显示隐藏旧的注释

khaoula Oueslati 2022年4月19日

你好 @Torsten ，我不知道如何实现使用Lsqcurvefit，因为我不太熟悉matlab，我正在努力使用它来解决我的问题。我尝试使用fminunc，但错误不断出现，我不知道如何解决他们。

这是我使用fminunc的代码

                                   函数dDP = Mymodel (t, DP)
                                  
                                   %已知参数
                                  
                                   E = 111年e3;
                                  
                                   R = 8.314;
                                  
                                   T = 371;
                                  
                                   dDP = - (A * (exp (e / R * T))) * (DP ^ 2);
                                  
                                   结束
                                  
                                   函数obj = objective(A)
                                  
                                   DP0 = 1000;
                                  
                                   %的时间点
                                  
                                   Ts =[1 2 3 4 5 6 7 8];
                                  
                                   (t, DP) =数值(@ (t, DP) Mymodel (t, DP), ts, DP0);
                                  
                                   DP_measured=[1000 700.32 580.42 408.20 317.38 281.18 198.15 100.12];
                                  
                                   = (DP-DP_measured) ^ 2;
                                  
                                   obj = (A)之和;
                                  
                                   结束
                                  
                                   A0 = 1 e8;
                                  
                                   乐趣= @目标;
                                  
                                   [A, fval] = fminunc(有趣,A0);
                                  
                                   disp ([”答:“num2str (A)])
                                  
                                   这些出现错误:
                                  
                                   odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
                                  
                                   错误在目标中(第5行)
                                  
                                   (t, DP) =数值(@ (t, DP) Mymodel (t, DP), ts, DP0);
                                  
                                   optimization_DP
                                  
                                   错误使用fminunc
                                  
                                   提供目标函数必须返回一个标量值。
                                  
                                   错误在optimization_DP(第3行)
                                  
                                   [A, fval] = fminunc(有趣,A0);

登录评论。

山姆翟 2022年4月14日

0
链接

直接链接到这个答案

https://ch.mathworks.com/matlabcentral/answers/1696395-solving-a-non-linear-ode-with-unknown-parameter#answer_942685

编辑:山姆翟 2022年4月14日

嗨 @khaoula Oueslati

这个控制方程已经给出了，你已经得到了

数据。

目标是想要找到一个．

从

数据，你可以估计

．接下来,

可以从微分方程：

现在,如果 R ， E 而且 T 都是已知的

可以从代数方程：

请核实一下。

如果

数据是均匀分布的，那么可以用这种方法进行估计

．

                              T = -pi:(2*pi/100):pi;
                             
                              X = sint;%假设Dp是正弦波
                             
                              Y = gradient(x)/(2*pi/100);%估计dotDp，余弦波预期
                             
                              情节(t x,“线宽”， 1.5 t, y，“线宽”, 1.5)
                             
                              网格在
                             
                              包含(“t”）
                             
                              ylabel (x(t)和x " (t)）
                             
                              传奇('x(t) = sint '，x " (t) = cos(t)，“位置”，“西北”）