[已解决] 利用图像拟合杨拉普拉斯微分方程组中的未知参数

[复制链接]
雨声之中 发表于 2021-9-17 16:43:09
10 财富积分
本帖最后由 雨声之中 于 2021-9-20 15:30 编辑

由于课题研究需要,临时接触matlab编程问题,希望通过图片形状得到杨拉普拉斯方程中的某一个参数,已经完成的部分是通过图片将图像坐标点提取出来了,但后续的参数拟合过程没有找到合适的参考代码所以迟迟没有进展,求有能力的朋友多多帮忙,修改后的代码在楼下更新,但仍有问题
杨-拉普拉斯微分方程组
dφ/ds=-sinφ/x+2/R-(998-1.29)*9.8*z/k(1)
dx/ds=cosφ
dz/ds=sinφ
0=x(s=0)=z(s=0)=φ(s=0)
R=0.004,当s=0时,(1)式变成dφ/ds=1/R
希望能通过图片得出k的值
修改别人的代码如下,但没有办法进行计算,
报错
错误使用 fmincon (第 232 行)
Invalid datatype. Options argument must be created with OPTIMOPTIONS.
出错 KineticsEst5 (第 43 行)
[k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],lb,ub,[],[],x0,yexp);
我认为主要还是拟合程序有问题,希望各位能提出宝贵的意见和建议,谢谢
  1. function KineticsEst5
  2. % 动力学ODE方程模型的参数估计
  3. %
  4. %
  5. % [Ref]:Frerich Keil, et al. ed., Scientific computing in chemical
  6. %   engineering II,1999 (P.351)
  7. %
  8. %   The variables  y here are y(1)=x1, y(2)=x4, y(3)=x5,y(4)=x6 .
  9. close all;clear all;clc;
  10. % 提取图像轮廓,提取图像边缘
  11. I= imread('333333.png');
  12. c = im2bw(I,graythresh(I));
  13. %figure;
  14. %subplot(131);
  15. b = edge(c,'canny');
  16. [u,v] = find(b);   %返回边界矩阵b中非零元素的位置
  17. xp = v;  %行值v赋给xp
  18. yp = u;  %列值u赋给yp
  19. x0 = mean([min(xp),max(xp)]);  %x0为行值的均值
  20. y0 = mean([min(yp),max(yp)]);  %y0为列值得均值
  21. xp1 = xp-x0;
  22. yp1 = yp-y0;
  23. [cita,r] = cart2pol(xp1,yp1);
  24. q = sortrows([cita,r]);  %从r列开始比较数值并按升序排序
  25. cita = q(:,1);  %赋角度值
  26. r = q(:,2);  %赋半径模值
  27. %subplot(132);polar(cita,r);  %画极坐标下的轮廓图
  28. [x,y] = pol2cart(cita,r);
  29. c=min(y+y0)
  30. y2 = x*0.0082317;%乘以像素点代表长度
  31. y1 = (y+y0-c)*0.0082317;%乘以像素点代表长度
  32. subplot(133);plot(y2,y1);axis equal;
  33. k0 = [0.0724];         % 参数初值
  34. lb = [0];                   % 参数下限
  35. ub = [+inf];    % 参数上限
  36. x0 = [0 0 0 0];
  37. ExpData=zeros(3332,4);
  38. ExpData(:,3)=y2;
  39. ExpData(:,4)=y1;
  40. %KineticsData1;
  41. yexp=ExpData;                  % yexp: 实验数据[x1        x4        x5        x6]
  42. % 使用函数fmincon()进行参数估计
  43. [k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],lb,ub,[],[],x0,yexp);
  44. fprintf('n使用函数fmincon()估计得到的参数值为:n')
  45. fprintf('tk1 = %.4fn',k(1))
  46. fprintf(' The sum of the squares is: %.1enn',fval)
  47. k_fmincon = k;
  48. %使用函数lsqnonlin()进行参数估计
  49. [k,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);      
  50. ci = nlparci(k,residual,jacobian);
  51. fprintf('nn使用函数lsqnonlin()估计得到的参数值为:n')
  52. Output
  53. % 以函数fmincon()估计得到的结果为初值,使用函数lsqnonlin()进行参数估计
  54. k0 = k_fmincon;
  55. [k,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);      
  56. ci = nlparci(k,residual,jacobian);
  57. fprintf('nn以fmincon()的结果为初值,使用函数lsqnonlin()估计得到的参数值为:n')
  58. Output
  59. % ------------------------------------------------------------------
  60. function f = ObjFunc4Fmincon(k,x0,yexp)

  61. tspan = [0.00 : 0.01 : 0.20];

  62. [t x] = ode45(@KineticEqs,tspan,x0,[],k);

  63. y(:,1) = x(:,1);

  64. y(:,2:4) = x(:,4:6);

  65. f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2)+ sum((y(:,3)-yexp(:,3)).^2) + sum((y(:,4)-yexp(:,4)).^2);
  66. % ------------------------------------------------------------------
  67. % ------------------------------------------------------------------
  68. function dxdt = KineticEqs(t,x,k)
  69. k=k(1);
  70. if t==0
  71.     dxdt=[250 0 0]%250为1/r顶点曲率
  72. else
  73.     P=(1e-6)/3*150*9.8*(-cos(x(1)))*((2-cos(pi*90/180))*(1+cos(pi*90/180))^2*(2500-998)+(2+cos(pi*90/180))*(1-cos(pi*90/180))^2*2500)
  74.     ;%修改接触角、覆盖率
  75.     dxdt =[
  76.     (-sin(x(1))/x(2)+2/0.00128-(998-1.29)*9.8*x(3)/k(1)+P/k(1))
  77.     (cos(x(1)))
  78.     (sin(x(1)))
  79.     ]
  80. end
复制代码


333333.png

最佳答案


TouAkira 发表于 2021-9-19 10:52:22
问题太多了…
最基本的,既然你的微分方程组中
dx/ds=cosφ
dz/ds=sinφ
那么,怎么能在KineticEqs中设置 t == 0 时,dxdt=[250 0 0] 呢? 这不符合 cos φ² + sin φ² == 1的基本三角函数公式。

还有,微分方程组里面就三个方程,你却在非线性拟合的目标函数中调用
[t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,2:4) = x(:,4:6);
x算出来就三列,不可能有x(:,4:6)的
回复此楼

4 条回复


TouAkira 发表于 2021-9-19 10:52:22
问题太多了…
最基本的,既然你的微分方程组中
dx/ds=cosφ
dz/ds=sinφ
那么,怎么能在KineticEqs中设置 t == 0 时,dxdt=[250 0 0] 呢? 这不符合 cos φ² + sin φ² == 1的基本三角函数公式。

还有,微分方程组里面就三个方程,你却在非线性拟合的目标函数中调用
[t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,2:4) = x(:,4:6);
x算出来就三列,不可能有x(:,4:6)的
回复此楼

雨声之中 发表于 2021-9-19 21:19:01
TouAkira 发表于 2021-9-19 10:52
问题太多了…
最基本的,既然你的微分方程组中
dx/ds=cosφ

您好,感谢您的回复,微分方程组应该是合理的,自变量为0的时候方程1是变成了等于1/R,也就是曲率,这个是方程自带的条件,如果没有的话因为全是0的初值是得不到结果的,但是代码确实是有问题的,因为之前有的地方我没仔细理解,现在我正在进一步修改,感谢您的回答,修改后如果有问题希望您能再指导一下,谢谢您:handshake

雨声之中 发表于 2021-9-20 15:27:30
本帖最后由 雨声之中 于 2021-9-20 17:27 编辑

修改后的代码参考论坛大佬TouAkira 的代码,分成两部分,一部分为提取轮廓数据并输出,一部分为拟合,提取轮廓如下
  1. clear; clc; close all;
  2. % 提取图像轮廓,提取图像边缘
  3. I= imread('3333331.png');
  4. c = im2bw(I,graythresh(I));
  5. b = edge(c,'canny');
  6. [u,v] = find(b);   %返回边界矩阵b中非零元素的位置
  7. xp = v;  %行值v赋给xp
  8. yp = u;  %列值u赋给yp
  9. x0 = mean([min(xp),max(xp)]);  %x0为行值的均值
  10. y0 = mean([min(yp),max(yp)]);  %y0为列值得均值
  11. xp1 = xp-x0;
  12. yp1 = yp-y0;
  13. [cita,r] = cart2pol(xp1,yp1);
  14. q = sortrows([cita,r]);  %从r列开始比较数值并按升序排序
  15. cita = q(:,1);  %赋角度值
  16. r = q(:,2);  %赋半径模值
  17. %subplot(132);polar(cita,r);  %画极坐标下的轮廓图
  18. [x,y]= pol2cart(cita,r);
  19. e=min(y+y0)
  20. x = x*0.0082317;%乘以像素点代表长度
  21. y = (y+y0-e)*0.0082317;%乘以像素点代表长度
  22. d=[x,y];
  23. [m,n]=find(x>=0);
  24. y2=d(m,1);
  25. y1=d(m,2);
  26. tspan=zeros(m(end)-m(1)+1,1);
  27. for i=2:length(tspan)
  28. tspan(i)=tspan(i-1)+sqrt((y2(i)-y2(i-1))^2+(y1(i)-y1(i-1))^2);
  29. end
  30. c1exp=NaN*ones(m(end)-m(1)+1,1);
  31. c2exp = y2; %实验数据
  32. c3exp = y1;
  33. data=[tspan,c1exp, c2exp,c3exp];
  34. save data2.txt -ascii data;

复制代码

拟合代码
  1. function nihe
  2. clear; clc; close all;
  3. I=load ('data2.txt')
  4. x=I(:,3);
  5. y=I(:,4);
  6. subplot(133);plot(x,y);axis equal;
  7. k0 = 0.0724 * ones( 1, 1 );
  8. lb = zeros( 1, 1 );
  9. ub = inf * ones( 1, 1);
  10. cq0 = [ 0,0,0 ];
  11. tspan=I(:,1);
  12. c1exp=I(:,2);
  13. c2exp = x; %实验数据
  14. c3exp = y;
  15. cexp = [ c1exp, c2exp,c3exp ];
  16. % 使用函数fmincon()进行参数估计
  17. MyFMC_Opt = optimoptions( @fmincon, 'Display', 'off', 'FinDiffType', 'central' );
  18. [ k, fval ] = fmincon( @( K ) ObjFun4Fmincon( K, cq0, cexp ), k0, [], [], [], [], lb, ub, [], MyFMC_Opt );
  19. fprintf( [ '\tk = ', repmat('\t%.4f\t', [ 1, numel( k ) ] ), '\n' ], k );
  20. fprintf( 'The sum of the squares is:%.1e\n\n', fval );
  21. k_fmincon = k;
  22. % 以函数fmincon()估计得到的结果为初值,使用函数lsqnonlin()进行参数估计
  23. k0 = k_fmincon;
  24. [ k, ~, residual, ~, ~, ~, jacobian ] = lsqnonlin( @( K ) ObjFunc4LNL( K, cq0, cexp ), k0, lb, ub, [] );
  25. ci = nlparci( k, residual, jacobian );
  26. fprintf( '\n\n以fmincon()的结果为初值,使用函数lsqnonlin()估计得到的参数值为:\n' );
  27. arrayfun( @( ii ) fprintf( [ '\tk', num2str( ii ), ' = %.4f, 95%% Confidence Interval ( %.4f, %.4f )\n' ], k( ii ), ci( ii, 1 ), ci( ii, 2 ) ,[ 1: 1 : numel( k ) ] );
  28. [ t, c ] = ode45( @KineticEqs, tspan(1) : 1 : tspan(end), cq0, [], k_fmincon );
  29. plot(c2exp,c3exp,'ko',c(:,2),c(:,3),'k--')
  30. end
  31. %---------------------------------------------------------------------
  32. function f = ObjFun4Fmincon( k, cq0, cexp )
  33. I=load ('data2.txt')
  34. tspan=I(:,1);
  35. [ ~, c ] = ode45( @KineticEqs, tspan, cq0, [], k );
  36. Index_1 = find( ~isnan( cexp( :, 1 ) ) );
  37. Index_2 = find( ~isnan( cexp( :, 2 ) ) );
  38. Index_3 = find( ~isnan( cexp( :, 3 ) ) );
  39. f = sum( ( c( Index_1, 1 ) - cexp( Index_1, 1 ) ).^2 ) +sum( ( c( Index_2, 2 ) - cexp( Index_2, 2 ) ).^2)+sum( ( c( Index_3, 3 ) - cexp( Index_3, 3 ) ).^2 );
  40. end
  41. %---------------------------------------------------------------------
  42. function f = ObjFunc4LNL( k, cq0, cexp )
  43. I=load ('data2.txt')
  44. tspan=I(:,1);
  45. [ ~, c ] = ode45( @KineticEqs, tspan, cq0, [], k );
  46. Index_1 = find( ~isnan( cexp( :, 1 ) ) );
  47. Index_2 = find( ~isnan( cexp( :, 2 ) ) );
  48. Index_3 = find( ~isnan( cexp( :, 3 ) ) );
  49. f1 = c( Index_1, 1 ) - cexp( Index_1, 1 );
  50. f2 = c( Index_2, 2 ) - cexp( Index_2, 2 );
  51. f3 = c( Index_3, 3 ) - cexp( Index_3, 3 );
  52. f = [ f1; f2;f3 ];
  53. end

  54. %---------------------------------------------------------------------
  55. function dcdt = KineticEqs( t, c, k )
  56. k=k(1);
  57. if t==0
  58.     dcdt=[686.34;);(cos(c(1)));(sin(c(1)))]%686.34为1/r顶点曲率
  59. else
  60.     %P=(1e-6)/3*150*9.8*(-cos(c(1)))*((2-cos(pi*90/180))*(1+cos(pi*90/180))^2*(2500-998)+(2+cos(pi*90/180))*(1-cos(pi*90/180))^2*2500);
  61. %修改接触角、覆盖率
  62.     dcdt =[(-sin(c(1))/c(2)+2/0.001457-(998-1.29)*9.8*c(3)/k(1);(cos(c(1)));(sin(c(1)))]
  63. end
  64. end
复制代码

但工作台一直在输出dcdx的值,并没有输出k的值,而且矩阵c的值明显错误,希望有能力的大佬们能帮忙指导一下,感谢
3333331.png

雨声之中 发表于 2021-9-20 16:45:24
TouAkira 发表于 2021-9-19 10:52
问题太多了…
最基本的,既然你的微分方程组中
dx/ds=cosφ

您好,我参照您回复其他人的帖子中的代码另外修改了一下,但最终并没有k的值输出,如果您有时间能麻烦帮忙看一下吗,非常感谢
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表