• MATLAB：拟合与插值

     一、关于多项式的基本操作

    若要求非线性方程的根，则采用fzero, fminbnd函数  

     

    二、多项式拟合

    1. clc, clear
    2. x=0:0.2:10; y=0.25*x+20*sin(x);
    3. plot(x,y,'k.','MarkerSize',15)
    4. grid on; 
    5. hold on
    6. [p1,s1,mu1]=polyfit(x,y,3); %3阶多项式拟合
    7. y1=polyval(p1,x,s1,mu1);
    8. [p2,s2,mu2]=polyfit(x,y,5); %5阶多项式拟合
    9. y2=polyval(p2,x,s2,mu2);
    10. [p3,s3,mu3]=polyfit(x,y,8); %8阶多项式拟合
    11. y3=polyval(p3,x,s3,mu3);
    12. plot(x,y1,'c-', x,y2, 'r-', x,y3, 'b-');
    13. xlabel('X'); 
    14. ylabel('Y');
    15. legend('原始数据','3阶多项式拟合', '5阶多项式拟合', '8阶多项式拟合','Location','best');
    16. legend('boxoff')
    17. title('不同次数拟合曲线对比图')

    多项式评价和置信区间的估计：

    1. clc,clear
    2. x=0:0.2:10; 
    3. y=0.25*x+20*sin(x);
    4. [p,s,mu]=polyfit(x,y,6); %6阶多项式拟合
    5. [Y,DELTA] = polyconf(p,x,s,0.05,mu);
    6. fill([x,fliplr(x)],[Y-DELTA,fliplr(Y+DELTA)],[0.95,0.899,0.9230]) %RGB
    7. hold on 
    8. plot(x,y,'k.','MarkerSize',15)
    9. plot(x,y,'b-')
    10. plot(x,Y-DELTA,'m--')
    11. plot(x,Y+DELTA,'m--')
    12. grid on 
    13. title('6次多项式拟合及置信区间')
    14. xlabel('X')
    15. ylabel('Y')

    曲线拟合案例分析： 

    1. clc,clear
    2. [Data,Textdata] = xlsread('matlab视频数据\examp41.xls',1,'B3:C61');
    3. x = 1:59;
    4. plot(x,Data,'r.','MarkerSize',15)
    5. % xtick坐标刻度
    6. % numel数组元素的个数（或用length）
    7. % xticklabel坐标显示的字符串
    8. set(gca,'XTick',1:8:numel(x),'xticklabel',Textdata(1:8:end))
    9. xlabel('时间')
    10. ylabel('食品零售价格分类指数')
    11. title('我国2007年1月至2011年11月食品零售价格分类指数数据')
    12. x = x';
    13.  
    14. hold on 
    15. [p4,s4,mu4] = polyfit(x,Data,4);
    16. y4 = polyval(p4,x,s4,mu4);
    17. plot(x,y4,'k:','LineWidth',1.5)
    18.  
    19. [p6,s6,mu6] = polyfit(x,Data,6);
    20. y6 = polyval(p6,x,s6,mu6);
    21. plot(x,y4,'m--','LineWidth',1.5)
    22.  
    23. [p9,s9,mu9] = polyfit(x,Data,9);
    24. y9 = polyval(p9,x,s9,mu9);
    25. plot(x,y9,'b.-','LineWidth',1.5)
    26.  
    27. [p11,s11,mu11] = polyfit(x,Data,11);
    28. y11 = polyval(p11,x,s11,mu11);
    29. plot(x,y11,'g-','LineWidth',1.5)
    30.  
    31. legend('原始数据','4次拟合','6次拟合','9次拟合','11次拟合'); 
    32. legend('boxoff')

    自定义函数拟合：

    1. function y = nlinfun(beta,x)
    2.     a = beta(1);
    3.     b = beta(2);
    4.     y = a + (0.49 - a).*exp(-b*(x-8));
    5. end
    1. clc,clear,close
    2. data = xlsread('matlab视频数据\nlinfit_data.xlsx');
    3. x = data(1,:)';
    4. y = data(2,:)';
    5.  
    6. beta0 = [1,1];
    7. beta = nlinfit(x,y,'nlinfun',beta0);
    8. yp = nlinfun(beta,x);
    9. plot(x,y,'r.','MarkerSize',15)
    10. hold on 
    11. grid on
    12. xlabel('时间T')
    13. ylabel('氯气积分Y')
    14. title('化工生产中氯气积分随时间下降拟合曲线')
    15. plot(x,yp,'b-','LineWidth',2)
    16.

    三、一维数据插值 

     一维插值函数插值方法对比：

    1. clc,clear
    2. x = 0:10; 
    3. y = sin(x); 
    4. xi = 0:0.1:10; %xi表示插值点
    5. strmod = {'nearest','linear','spline','cubic'}; % 将插值方法定义为单元数组
    6. strlb = {'(a) method = nearest', '(b) method = linear','(c) method = spline', '(d) method = cubic'}; % 将X轴标识为单元数组
    7. for i = 1:4
    8.     yi = interp1(x,y,xi,strmod{i}); %一维插值
    9.     subplot(2,2,i); %生成子图
    10.     plot(x,y,'ro','MarkerFaceColor','r');
    11.     hold on
    12.     grid on
    13.     plot(xi,yi,'b--','LineWidth',1.5)
    14.     title(strlb(i)) %对每个字图添加标题
    15. end

    案例：环境温度数据插值

    1. clc,clear
    2. x = 0:2:24;
    3. y = [12 9 9 10 18 24 28 27 25 20 18 15 13];
    4. xi = 0:24/1440:24;
    5.  
    6. yisp = spline(x,y,xi);
    7. % P = spline(x,y);
    8. % yisp = interp1(x,y,xi,'spline');
    9. subplot(2,2,1)
    10. plot(x,y,'bo',xi,yisp,'r-')
    11. title('spline函数插值效果图')
    12. xlabel('24时间'); ylabel('随时间温度变化值')
    13. grid on
    14.  
    15. subplot(2,2,2)
    16. pcs = csape(x,y,'complete') %查看三次样条插值系数矩阵
    17. ycs = fnval(pcs,xi); %求插值
    18. plot(x,y,'bo',xi,ycs,'r-')
    19. title('csape函数插值效果图')
    20. xlabel('24时间'); ylabel('随时间温度变化值')
    21. grid on
    22.  
    23. subplot(2,2,3)
    24. %B样条插值， k为B样条阶次，一般选择4和5
    25. psp = spapi(4,x,y);
    26. yspa = fnval(psp,xi);
    27. plot(x,y,'bo',xi,yspa,'r-')
    28. grid on
    29. xlabel('24时间'); ylabel('随时间温度变化值')
    30. title('spapi函数插值效果图')
    31.  
    32. subplot(2,2,4)
    33. %三次光滑样条插值， p表示光滑程度，取值[0,1]
    34. ycsa = csaps(x,y,0.9,xi);
    35. plot(x,y,'bo',xi,ycsa,'r-')
    36. title('csaps函数插值效果图')
    37. xlabel('24时间'); ylabel('随时间温度变化值')
    38. grid on

    案例：轮船甲板面积

    1. clc,clear
    2. x = linspace(0,8.534,13);
    3. y = [0 0.914 5.060 7.772 8.717 9.083 9.144 9.083 8.992 8.687 7.376 2.073 0];
    4. x0 = 0:0.001:8.534; %插值点
    5. y1 = interp1(x,y,x0,'linear'); %线性插值
    6. y2 = interp1(x,y,x0,'spline'); %三次样条插值
    7. plot(x,y,'b.','Markersize',15); 
    8. hold on
    9. plot(x0,y1,'r--',x0,y2,'g-');
    10. S1 = trapz(y1)*0.001 %线性插值数值积分，计算梯形面积
    11. S2 = trapz(y2)*0.001 %三次样条插值数值积分

    外插估值：

    1. % cos(10.5)外插值计算， x的区间[0.10];利用函数interp1计算10.5的函数值。
    2. clc,clear
    3. x = 0:0.5:10;
    4. y = cos(x);
    5. x1 = 10.5;
    6. y1 = cos(x1);
    7. y2 = interp1(x,y,x1,'nearest','extrap');
    8. y3 = interp1(x,y,x1,'linear','extrap');
    9. y4 = interp1(x,y,x1,'spline','extrap');
    10. plot(x,y,x1,y1,'o',x1,y2,'>r',x1,y3,'b<',x1,y4,'r*');
    11. legend('solution','cos(10.5)','nearest','linear','spline');
    12. grid on
    13. title('各种外插方法对比图')
    14. legend('boxoff')

    四、高维数据插值 

    案例：山区地貌图

    1. % 已知某处山区地形选点测量坐标数据为：x、y和z。
    2. clc,clear
    3. x = 0:0.5:5; y = 0:0.5:6;
    4. z = [89 90 87 85 92 91 96 93 90 87 82;
    5.     92 96 98 99 95 91 89 86 84 82 84;
    6.     96 98 95 92 90 88 85 84 83 81 85;
    7.     80 81 82 89 95 96 93 92 89 86 86;
    8.     82 85 87 98 99 96 97 88 85 82 83;
    9.     82 85 89 94 95 93 92 91 86 84 88;
    10.     88 92 93 94 95 89 87 86 83 81 92;
    11.     92 96 97 98 96 93 95 84 82 81 84;
    12.     85 85 81 82 80 80 81 85 90 93 95;
    13.     84 86 81 98 99 98 97 96 95 84 87;
    14.     80 81 85 82 83 84 87 90 95 86 88;
    15.     80 82 81 84 85 86 83 82 81 80 82;
    16.     87 88 89 98 99 97 96 98 94 92 87];
    17. mesh(x,y,z) %绘制三维网格面
    18.  
    19. xi = linspace(0,5,50);%加密横坐标数据到50个
    20. yi = linspace(0,6,60);%加密纵坐标数据到60个
    21. [xi1,yi1] = meshgrid(xi,yi);%生成网格数据
    22. [xi2,yi2] = ndgrid(xi,yi);
    23.  
    24. % ------- 调用interp2函数作三次样条插值-----------
    25. subplot(2,3,1);
    26. zi1 = interp2(x,y,z,xi1,yi1,'spline');
    27. mesh(xi1,yi1,zi1); 
    28. hold on
    29. [xx,yy] = meshgrid(x,y); %生成网格数据
    30. plot3(xx,yy,z+0.1,'bo') %原始数据用O绘制
    31. title('interp2函数绘图')
    32.  
    33. % ------- 调用csape函数作三次样条插值----------
    34. subplot(2,3,2);
    35. cp1 = csape({x,y},z');
    36. mesh(xi2,yi2,fnval(cp1,{xi,yi})); 
    37. hold on
    38. [xx,yy] = meshgrid(x,y); %生成网格数据
    39. plot3(xx,yy,z+0.1,'bo') %原始数据用O绘制
    40. title('csape函数绘图')
    41.  
    42. % ------调用csapi函数作三次样条插值---------
    43. subplot(2,3,3)
    44. zi2 = csapi({x,y},z',{xi,yi});
    45. mesh(xi2,yi2,zi2);
    46. hold on
    47. [xx,yy] = meshgrid(x,y); %生成网格数据
    48. plot3(xx,yy,z+0.1,'bo') %原始数据用O绘制
    49. title('csapi函数绘图')
    50.  
    51. % ------调用spapi函数作三次B样条插值--------
    52. subplot(2,3,4)
    53. sp1 = spapi({4,4},{x,y},z');
    54. mesh(xi2,yi2,fnval(sp1,{xi,yi}))
    55. hold on
    56. [xx,yy] = meshgrid(x,y); %生成网格数据
    57. plot3(xx,yy,z+0.1,'bo') %原始数据用O绘制
    58. title('spapi函数绘图')
    59.  
    60. subplot(2,3,5)
    61. zi3 = csaps({x,y},z',{0.2,0.9},{xi,yi});
    62. mesh(xi2,yi2,zi3);
    63. hold on
    64. [xx,yy] = meshgrid(x,y); %生成网格数据
    65. plot3(xx,yy,z+0.1,'bo') %原始数据用O绘制
    66. title('csaps函数绘图')
    67.  
    68. % ---------调用spaps函数作三次B样条插值-------------
    69. subplot(2,3,6)
    70. sp2 = spaps({x,y},z',{1e-3,0.01});
    71. mesh(xi2,yi2,fnval(sp2,{xi,yi}))
    72. hold on
    73. [xx,yy] = meshgrid(x,y); %生成网格数据
    74. plot3(xx,yy,z+0.1,'bo') %原始数据用O绘制
    75. title('spaps函数绘图')

    散乱点插值：

    案例：水道海底地貌图

    1. clc,clear
    2. x = [129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5];
    3. y = [7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5];
    4. z = -[4 8 6 8 6 8 8 9 9 8 8 9 4 9];
    5. [cx,cy] = meshgrid(75:5:200,-90:5:150);
    6. cz = griddata(x,y,z,cx,cy,'cubic');
    7. figure(1)
    8. mesh(cx,cy,cz); %绘制三维网格图
    9. title('某水道的海底地貌图')
    10. view(-60,30);
    11. figure(2)
    12. [c,h] = contourf(cx,cy,cz,[-5,-5],'k'); %绘制等高线
    13. set(h,'ShowText','on')
    14. grid on
    15. title('船只禁入区域图')

    案例：城区土壤地质环境调查

    1. clc,clear
    2. data = xlsread('matlab视频数据\cumcm2011A.xls',1,'B4:D322');
    3. x = data(:,1);
    4. y = data(:,2);
    5. z = data(:,3);
    6. cd = xlsread('matlab视频数据\cumcm2011A.xls',2,'C4:C322');
    7. xd = linspace(min(x),max(x),60);
    8. yd = linspace(min(y),max(y),60);
    9. [xi,yi] = meshgrid(xd,yd);
    10.  
    11. % ------调用griddata函数作散乱节点插值-------
    12. zi = griddata(x,y,z,xi,yi);
    13. cdi = griddata(x,y,cd,xi,yi);
    14. surf(xi,yi,zi,cdi)  % 前三维度绘制空间曲面，第四维度用颜色表示
    15. shading interp
    16. xlabel('X'); 
    17. ylabel('Y'); 
    18. zlabel('Z（griddata）');
    19. colorbar
    20.  
    21. % ------------调用TriScatteredInterp函数作散乱节点插值------------
    22. F = TriScatteredInterp(x,y,z);
    23. zi2 = F(xi,yi);
    24. Fcd = TriScatteredInterp(x,y,cd);
    25. cdi2 = Fcd(xi,yi);
    26. figure;
    27. surf(xi,yi,zi2,cdi2);
    28. shading interp;
    29. xlabel('X');
    30. ylabel('Y');
    31. zlabel('Z（TriScatteredInterp）');
    32. colorbar

     

     

  • 相关阅读:
    Milvus Cloud——LLM Agent 现阶段出现的问题
    DDR3 NATIVE接口
    2022牛客暑期多校训练营7(总结+补题)
    Python自学教程2：大牛们怎么写注释
    【操作系统】存储器管理：对换
    [acwing周赛复盘] 第 64 场周赛20220813
    C Primer Plus(6) 中文版 第6章 C控制语句：循环 6.5 for循环
    技术分享| gcc版本升级到5.2
    【DL】基于卷积神经网络的一维数据二分类预测
    scala中函数名 _(下划线)如何理解？
  • 原文地址：https://blog.csdn.net/2302_78896863/article/details/136747093