腐蚀电化学极化曲线拟合程序

曲线拟合及其MATLAB 程序

实验任务：用最小二乘法将一组离散数据进行曲线拟合，并用MATLAB 做出曲线。

实验要求：

第三次作业的最后一题求函数关系并作出曲线，此组x 和y 的实验数据如下表所示。 x i

1 2 3 4 5 6 7 8 9 10 y i 6.00 3.50 2.67 2.75 2.00 1.80 1.70 1.63 1.56 1.50

实验内容：

（1）在在MATLAB 工作窗口输入程序

>> x=[1 2 3 4 5 6 7 8 9 10];

y=[6.00 3.50 2.67 2.75 2.00 1.80 1.70 1.63 1.56 1.50];

plot(x,y,'r*'),

legend('实验数据(xi,yi)')

xlabel('x'), ylabel('y'),

title('实验数据点(xi,yi)的散点图')

（2）运行后屏幕显示数据的散点图（略）。

（3）编写下列MATLAB 程序计算)(x f 在),(i i y x 处的函数值，即输入程序 >> syms a1 a2 a3 a4;

x=[1 2 3 4 5 6 7 8 9 10];

fi=a1.*x.^3+ a2.*x.^2+ a3.*x+ a4

运行后屏幕显示关于a1, a2, a3和a4的线性方程组

fi =

[a1+a2+a3+a4, 8*a1+4*a2+2*a3+a4, 27*a1+9*a2+3*a3+a4, 64*a1+16*a2+4*a3+a4, 125*a1+25*a2+5*a3+a4,

216*a1+36*a2+6*a3+a4, 343*a1+49*a2+7*a3+a4,

512*a1+64*a2+8*a3+a4, 729*a1+81*a2+9*a3+a4,

1000*a1+100*a2+10*a3+a4]

编写构造误差平方和的MATLAB 程序

>> y=[6.00 3.50 2.67 2.75 2.00 1.80 1.70 1.63 1.56 1.50];

fi

=[a1+a2+a3+a4,8*a1+4*a2+2*a3+a4,27*a1+9*a2+3*a3+a4,64*a1+16*a2+4*a3+a4,125*a1+25*a2+5*a3+a4,216*a1+36*a2+6*a3+a4,343*a1+49*a2+7*a3+a4,512*a1+64*a2+8*a3+a4,729*a1+81*a2+9*a3+a4, 1000*a1+100*a2+10*a3+a4]; fy=fi-y; fy2=fy.^2;J=sum(fy.^2)

运行后屏幕显示误差平方和如下

J =

(a1+a2+a3+a4-6)^2+(8*a1+4*a2+2*a3+a4-7/2)^2+(27*a1+9*a2+3*a3+a4-267/100)^2+(64*a1+16*a2+4*a3+a4-11/4)^2+(125*a1+25*a2+5*a3+a4-2)^2+(216*a 1+36*a2+6*a3+a4-9/5)^2+(343*a1+49*a2+7*a3+a4-17/10)^2+(512*a1+64*a2+8*a3+a4-163/100)^2+(729*a1+81*a2+9*a3+a4-39/25)^2+(1000*a1+100*a2+10*a3+a4-3/2)^2

为求4321,,,a a a a 使J 达到最小，只需利用极值的必要条件0

=??k a J 。)4,3,2,1(=k ，得到关于4321,,,a a a a 的线性方程组，这可以由下面的MATLAB 程序完成，即输入程序

>> syms a1 a2 a3 a4

J=(a1+a2+a3+a4-6)^2+(8*a1+4*a2+2*a3+a4-7/2)^2+(27*a1+9*a2+3*a3+a4

-267/100)^2+(64*a1+16*a2+4*a3+a4-11/4)^2+(125*a1+25*a2+5*a3+a4-2)^2+( 216*a1+36*a2+6*a3+a4-9/5)^2+(343*a1+49*a2+7*a3+a4-17/10)^2+(512*a1+6 4*a2+8*a3+a4-163/100)^2+(729*a1+81*a2+9*a3+a4-39/25)^2+(1000*a1+100* a2+10*a3+a4-3/2)^2;

Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3); Ja4=diff(J,a4);

Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3), Ja41=simple(Ja4),

运行后屏幕显示J分别对a1, a2 ,a3 ,a4的偏导数如下

Ja11 =6050*a4+3956810*a1+441650*a2+50666*a3-497579/50

Ja21 =770*a4+441650*a1+50666*a2+6050*a3-66681/50

Ja31 =110*a4+50666*a1+6050*a2+770*a3-10679/50

Ja41 = 20*a4+6050*a1+770*a2+110*a3-2511/50

解线性方程组Ja11 =0，Ja21 =0，Ja31 =0，Ja41 =0，输入下列程序

>>A=[3956810, 441650, 50666, 6050; 441650, 50666, 6050, 770; 50666, 6050, 770, 110; 6050, 770, 110, 20];

B=[497579/50, 66681/50, 10679/50, 2511/50];

C=B/A, f=poly2sym(C)

运行后屏幕显示拟合函数f及其系数C如下

C = -0.0176 0.3744 -2.6452 7.9627

f=-5066427105276917/288230376151711744*x^3+6744376644764349/1801439 8509481984*x^2-2978180218619021/1125899906842624*x+112064********* 3/140737488355328

故所求的拟合曲线为

f(x)= -0.0176 x3 + 0.3744 x2 -2.6452 x+ 7.9627 .

（4）编写下面的MATLAB程序估计其误差，并作出拟合曲线和数据的图形，输入程序

>> xi =[1 2 3 4 5 6 7 8 9 10];

y=[6.00 3.50 2.67 2.75 2.00 1.80 1.70 1.63 1.56 1.50];

n=length(xi);

f= -0.0176.*xi.^3+0.3744.*xi.^2-2.6452.*xi +7.9627;

x=1: 1: 10;

F= -0.0176.*xi.^3+0.3744.*xi.^2 -2.6452.*xi +7.9627; fy=abs(f-y); fy2=fy.^2; Ew=max(fy),

E1=sum(fy)/n, E2=sqrt((sum(fy2))/n)

plot(xi,y,'r*'), hold on, plot(x,F,'b-'), hold off legend('数据点(xi,yi)','拟合曲线y=f(x)'), xlabel('x'), ylabel('y'),

title('实验数据点(xi,yi)和拟合曲线y=f(x)的图形') 运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差Ew ，平均误差E1和均方根误差E2及其数据点),(i i y x 和拟合曲线y=f(x)的图形如下所示。 Ew = 0.5291

E1 = 0.2163

E2 = 0.2763