数值分析上机作业
数值分析上机作业(总)
数值分析上机实验一、解线性方程组直接法（教材49页14题）追赶法程序如下：function x=followup(A,b)n = rank(A);for(i=1:n)if(A(i,i)==0)disp('Error: 对角有元素为0');return;endend;d = ones(n,1);a = ones(n-1,1);c = ones(n-1);for(i=1:n-1)a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i);endd(n,1) = A(n,n);for(i=2:n)d(i,1)=d(i,1) - (a(i-1,1)/d(i-1,1))*c(i-1,1);b(i,1)=b(i,1) - (a(i-1,1)/d(i-1,1))*b(i-1,1);endx(n,1) = b(n,1)/d(n,1);for(i=(n-1):-1:1)x(i,1) = (b(i,1)-c(i,1)*x(i+1,1))/d(i,1);end主程序如下：function zhunganfaA=[2 -2 0 0 0 0 0 0;-2 5 -2 0 0 0 0 0;0 -2 5 -2 0 0 0 0;0 0 -2 5 -2 0 0 0;0 0 0 -2 5 -2 0 0;0 0 0 0 -2 5 -2 0;0 0 0 0 0 -2 5 -2;0 0 0 0 0 0 -2 5];b=[220/27;0;0;0;0;0;0;0];x=followup(A,b)计算结果：x =8.14784.07372.03651.01750.50730.25060.11940.0477二、解线性方程组直接法（教材49页15题）程序如下：function tiaojianshu(n)A=zeros(n);for j=1:1:nfor i=1:1:nA(i,j)=(1+0.1*i)^(j-1);endendc=cond(A)d=rcond(A)当n=5时c =5.3615e+005d =9.4327e-007当n=10时c =8.6823e+011d =5.0894e-013当n=20时c =3.4205e+022d =8.1226e-024备注：对于病态矩阵A来说，d为接近0的数；对于非病态矩阵A来说，d为接近1的数。
数值分析上机作业(MATLAB)
将系数矩阵 A 分解为：A=L+U+D
Ax=b
⇔ (D + L +U)x = b ⇔ Dx = −(L + U )x + b ⇔ x = −D −1(L + U )x + D −1b x(k +1) = −D −1 (L + U ) x(k ) + D −1b
输入 A，b 和初始向量 x
迭代矩阵 BJ , BG
否
ρ(B) < 1?
按雅各比方法进行迭代
否
|| x (k+1) − x(k) ||< ε ?
按高斯-塞德尔法进行迭代
否
|| x(k+1) − x (k ) ||< ε ?
输出迭代结果
图 1 雅各布和高斯-赛德尔算法程序流程图
1.2 问题求解
按图 1 所示的程序流程，用 MATLAB 编写程序代码，具体见附录 1。解上述三个问题 如下
16
-0.72723528355328
0.80813484897616
0.25249261987171
17
-0.72729617968010
0.80805513082418
0.25253982509100
18
-0.72726173942623
0.80809395746552
0.25251408253388
0.80756312717373
8
-0.72715363032573
0.80789064377799
9
-0.72718652854079
数值分析大作业
数值分析上机作业（一）一、算法的设计方案1、幂法求解λ1、λ501幂法主要用于计算矩阵的按模最大的特征值和相应的特征向量，即对于|λ1|≥|λ2|≥.....≥|λn|可以采用幂法直接求出λ1，但在本题中λ1≤λ2≤……≤λ501，我们无法判断按模最大的特征值。
但是由矩阵A的特征值条件可知|λ1|和|λ501|之间必然有一个是最大的，通过对矩阵A使用幂法迭代一定次数后得到满足精度ε=10−12的特征值λ0，然后在对矩阵A做如下的平移：B=A-λ0I由线性代数(A-PI)x=(λ-p)x可得矩阵B的特征值为：λ1-λ0、λ2-λ0…….λ501-λ0。
对B矩阵采用幂法求出B矩阵按模最大的特征值为λ∗=λ501-λ0，所以λ501=λ∗+λ0，比较λ0与λ501的大小，若λ0>λ501则λ1=λ501，λ501=λ0；若λ0<λ501，则令t=λ501，λ1=λ0，λ501=t。
求矩阵M按模最大的特征值λ的具体算法如下：任取非零向量u0∈R nηk−1=u T(k−1)∗u k−1y k−1=u k−1ηk−1u k=Ay k−1βk=y Tk−1u k(k=1,2,3……)当|βk−βk−1||βk|≤ε=10−12时，迭终终止，并且令λ1=βk2、反幂法计算λs和λik由已知条件可知λs是矩阵A 按模最小的特征值，可以应用反幂法直接求解出λs。
使用带偏移量的反幂法求解λik，其中偏移量为μk=λ1+kλ501−λ140(k=1,2,3…39)，构造矩阵C=A-μk I，矩阵C的特征值为λik−μk，对矩阵C使用反幂法求得按模最小特征值λ0，则有λik=1λ0+μk。
求解矩阵M按模最小特征值的具体算法如下：任取非零向量u 0∈R n ηk−1= u T (k−1)∗u k−1y k−1=u k−1ηk−1 Au k =y k−1βk =y T k−1u k (k=1,2,3……)在反幂法中每一次迭代都要求解线性方程组Au k =y k−1，当K 足够大时，取λn =1βk 。
数值分析上机题目
数值分析上机题目4(总21页) --本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--实验一实验项目：共轭梯度法求解对称正定的线性方程组 实验内容：用共轭梯度法求解下面方程组(1) 123421003131020141100155x x x x -⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪=⎪ ⎪ ⎪-- ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ 迭代20次或满足()(1)1110k k x x --∞-<时停止计算。
编制程序：储存m 文件function [x,k]=CGmethod(A,b)n=length(A);x=2*ones(n,1);r=b-A*x;rho=r'*r; k=0;while rho>10^(-11) & k<1000 k=k+1; if k==1 p=r; elsebeta=rho/rho1; p=r+beta*p; end w=A*p;alpha=rho/(p'*w); x=x+alpha*p; r=r-alpha*w; rho1=rho;rho=r'*r; end运行程序： clear,clcA=[2 -1 0 0;-1 3 -1 0;0 -1 4 -1;0 0 -1 5]; b=[3 -2 1 5]'; [x,k]=CGmethod(A,b)运行结果： x =(2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵， A[i,i]=4，A[i,i-1]=A[i-1,i]=-1，i=2,3,..,n b[1]=3, b[n]=3, b[i]=2，i=2,3,…,n-1迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。
编制程序：储存m 文件function [x,k]=CGmethod_1(A,b) n=length(A);x(1:n,1)=0;r=b-A*x;r1=r; k=0;while norm(r1,1)>10^(-7)&k<10^4 k=k+1; if k==1 p=r; elsebeta=(r1'*r1)/(r'*r);p=r1+beta*p; end r=r1; w=A*p;alpha=(r'*r)/(p'*w); x=x+alpha*p; r1=r-alpha*w; end运行程序： clear,clc n=1000; A=hilb(n); b=sum(A')';[x,k]=CGmethod(A,b)实验二1、 实验目的：用复化Simpson 方法、自适应复化梯形方法和Romberg 方法求数值积分。
数值分析报告上机题课后作业全部-东南大学
实用标准文案文档大全上机作业题报告2015.1.9 USER1.Chapter 11.1题目设S N =∑1j 2−1N j=2，其精确值为)11123(21+--N N 。
（1）编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算S N 的通用程序。
（2）编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算S N 的通用程序。
（3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。
（编制程序时用单精度） （4）通过本次上机题，你明白了什么？1.2程序1.3运行结果1.4结果分析按从大到小的顺序，有效位数分别为：6，4，3。
按从小到大的顺序，有效位数分别为：5，6，6。
可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。
当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。
因此，采取从小到大的顺序累加得到的结果更加精确。
2.Chapter 22.1题目（1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。
（2）给定方程03)(3=-=x xx f ,易知其有三个根3,0,3321=*=*-=*x x x○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。
试确定尽可能大的δ。
○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。
（3）通过本上机题，你明白了什么？2.2程序2.3运行结果（1）寻找最大的δ值。
算法为：将初值x0在从0开始不断累加搜索精度eps，带入Newton迭代公式，直到求得的根不再收敛于0为止，此时的x0值即为最大的sigma值。
运行Find.m，得到在不同的搜索精度下的最大sigma值。
数值分析上机作业
第一章第二题(1) 截断误差为104-时：k=1;n=0;m=0;x=0;e=1e-4;while k==1x1=x+(-1)^n/(2*n+1);if abs(x-x1)<ey=4*x1;m=n+1;break;endx=x1;k=1;n=n+1;endformat longy,my =3.141792613595791m =5001（2）截断误差为108-时：k=1;n=0;m=0;x=0;e=1e-8;while k==1x1=x+(-1)^n/(2*n+1);if abs(x-x1)<ey=4*x1;m=n+1;break;endx=x1;k=1;n=n+1;endformat longy,my =3.141592673590250m =50000001由以上计算可知，截断误差小于104-时，应取5001项求和，π=3.141792613595791；截断误差小于108-时，应取50000001项求和，π=3.141592673590250。
第二章第二题a=[0 -2 -2 -2 -2 -2 -2 -2];b=[2 5 5 5 5 5 5 5];c=[-2 -2 -2 -2 -2 -2 -2 0];v=220;r=27;d=[v/r 0 0 0 0 0 0 0];n=8;for i=2:na(i)=a(i)/b(i-1);b(i)=b(i)-c(n-1)*a(i);d(i)=d(i)-a(i)*d(i-1);end;d(n)=d(n)/b(n);for i=n-1:-1:1d(i)=(d(i)-c(i)*d(i+1));end;I=d'I =1.0e+002 *1.490717294184090.704617906351300.311568212434910.128623612390290.049496991380330.017168822994210.004772412363470.00047741598598第三章第一题(1)Jacobi迭代法：b=[12;-27;14;-17;12]x = [0;0;0;0;0;]k = 0;r = 1;e=0.000001A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;] D = diag(diag(A));B = inv(D)*(D-A);f = inv(D)*b;p = max(abs(eig(B)));if p >= 1'迭代法不收敛'returnendwhile r >ex0 = x;x = B*x0 + f;k = k + 1;r = norm (x-x0,inf);endxk计算结果：x =1.0000-2.00003.0000-2.00001.0000k =65(2) 高斯赛德尔迭代：A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;]x=[0;0;0;0;0];b=[12;-27;14;-17;12]c=0.000001L=-tril(A,-1)U=-triu(A,1)D=(diag(diag(A)))X=inv(D-L)*U*x+inv(D-L)*b;k=1;while norm(X-x,inf)>= cx=X;X=inv(D-L)*U*x+inv(D-L)*b;k=k+1;endXk计算结果：X =1.0000-2.00003.0000-2.00001.0000k =37(3) SOR：A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15] x=[0;0;0;0;0];b=[12;-27;14;-17;12]e=0.000001w=1.44;L=-tril(A,-1)U=-triu(A,1)D=(diag(diag(A)))X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*bn=1;while norm(X-x,inf)>=ex=X;X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*b;n=n+1;endXn计算结果：X =1.0000-2.00003.0000-2.00001.0000n =22由以上可知，共轭梯度法收敛速度明显快于Jacobi法和G-S法。
东南大学出版社第二版《数值分析》上机作业答案(前三章)
for (i=k+1;i<N;i++) // { lik=a[i][k]/a[k][k]; //实施消去过程，得到上三角系数增广矩阵 for (j=k;j<M;j++) // { a[i][j]=a[i][j]‐lik*a[k][j]; // } } } cout<<"经列主元高斯消去法得到的数组为:"<<endl; // for (b=0;b<N;b++) // { cout<<endl; //输出经过列主元消去法处理过的系数增广矩阵 for (c=0;c<M;c++) { cout<<setw(7)<<a[b][c]; // } } cout<<endl; double x[N]; // double s; int f,g; x[N‐1]=a[N‐1][M‐1]/a[N‐1][N‐1]; // for (f=N‐2;f>=0;f‐‐) // { s=0; for (g=f+1;g<N;g++) //由上三角形的系数增广矩阵求出方程组的解 { s=s+a[f][g]*x[g]; // } x[f]=(a[f][N]‐s)/a[f][f]; // } cout<<"方程组的解为:"<<endl; for (b=0;b<N;b++) //输出方程组的解 {
1
当 n=10000 时，s3=0.7499 Press any key to continue （分析 S1 的 6 位数字中，有效位数为 4 位； S2 的所有数字都是有效数字。 ） 当 n=1000000 时，s1=‐14.2546 当 n=1000000 时，s2=‐14.2551 当 n=1000000 时，s3=0.749999 Press any key to continue （分析： S1 的 6 位数字中，没有有效数字； S2 的 6 位数字中，没有有效数字。 ） 由运行结果可知，当精度比较低时，按从大数开始累加到小数的计算结果的精度低于按从小数 累加到大数的计算结果的精度。 至于当 n=1000000 时，S1 和 S2 得出了负数结果，可能是由于循环次数过多，导致数据溢出， 从而得出错误结果。 习题 2 20．程序如下： //给定误差限为：0.5e‐6 //经过试算得当 delta 最大取道 0.7745966 时，迭代得到的根都收敛于 0 #include <iostream.h> #include <math.h> void main () { double x,u; int count=0; u=10.0; cout<<"请输入 x 的初值"<<endl; cin>>x; for (count=0;abs(u)>5;count++) { x=x‐(x*x*x‐3*x)/(3*(x*x‐1)); u=10000000*x; if(count>5000) { cout<<"迭代结果不收敛于 0！"<<endl; break; } } cout<<"x="<<x<<endl<<endl;
数值分析上机实验题参考
数值分析论文数值积分 一、问题提出选用复合梯形公式，复合Simpson 公式，Romberg 算法，计算I = dx x ⎰-4102sin 4 ()5343916.1≈II =dx x x ⎰1sin ()9460831.0,1)0(≈=I fI = dx xe x⎰+1024 I =()dx x x ⎰++1211ln 二、要求编制数值积分算法的程序；分别用两种算法计算同一个积分，并比较其结果；分别取不同步长()/ a b h -=n ，试比较计算结果（如n = 10, 20等）； ﹡给定精度要求ε，试用变步长算法，确定最佳步长﹡。
三、目的和意义深刻认识数值积分法的意义； 明确数值积分精度与步长的关系；根据定积分的计算方法，可以考虑二重积分的计算问题引言一、数值求积的基本思想实际问题当中常常需要计算积分，有些数值方法。
如微分方程和积分方程的求解，也都和积分计算相联系。
依据人们熟悉的微积分基本原理，对于积分I=⎰a b f(x)dx,只要找到被积函数f(x)和原函数F(x),便有下列牛顿-莱布尼茨公式：I=⎰a b f(x)dx=F(b)-F(a).但实际使用这种求积方法往往有困难，因为大量的被积函数，诸如x xsin，2xe-等，其原函数不能用初等函数表达，故不能用上述公式计算。
另外，当f(x)是由测量或数值计算给出的一张数据表时，牛顿-莱布尼茨公式也不能直接运用，因此有必要研究积分的数值计算问题。
二、数值积分代数精度数值求积方法是近似方法，为要保证精度，我们自然希望求积公式能对“尽可能多”的函数准确成立，就提出了所谓代数精度的概念。
如果某个求积公式对次数不超过m的多项式均能准确成立，但对m+1次多项式就不能准确成立，则称该求积公式具有m次代数精度。
三、复合求积公式为了提高精度，通常可以把积分区间分成若干子区间（通常是等分），再在每个子区间用低阶求积公式，即复化求积法，比如复化梯形公式与复化辛普森公式。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析上机实验报告选题：曲线拟合的最小二乘法指导老师：专业：学号：姓名：课题八 曲线拟合的最小二乘法一、问题提出从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。
在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。
二、要求1、用最小二乘法进行曲线拟合；2、近似解析表达式为()33221t a t a t a t ++=ϕ；3、打印出拟合函数()t ϕ，并打印出()j t ϕ与()j t y 的误差，12,,2,1 =j ；4、另外选取一个近似表达式，尝试拟合效果的比较；5、*绘制出曲线拟合图*。
三、目的和意义1、掌握曲线拟合的最小二乘法；2、最小二乘法亦可用于解超定线代数方程组；3、探索拟合函数的选择与拟合精度间的关系。
四、计算公式对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为∑==mj j j x a x y 0)()(ϕ特别的，取)(x j ϕ为多项式j j x x =)(ϕ (j=0, 1,…，m )则根据最小二乘法原理，可以构造泛函∑∑==-=ni mj i j j i m x a f a a a H 110))((),,,(ϕ令0=∂∂ka H(k=0, 1,…，m ) 则可以得到法方程⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ求该解方程组，则可以得到解m a a a ,,,10 ，因此可得到数据的最小二乘解∑=≈mj j j x a x f 0)()(ϕ曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。
曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。
五、结构程序设计在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。
5.1用一元三次多项式()33221t a t a t a t ++=ϕ进行拟合计算解析表达式系数： a1, a2, a3t=[0 5 10 15 20 25 30 35 40 45 50 55];y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64]; >> n=length(xi);f=0.34364.*10.^(-4)*x.^3-5.2156.*10.^(-3)*x.^2+0.26340.*x+0.017839;x=0:0.01:55;F=0.34364.*10.^(-4)*x.^3-5.2156.*10.^(-3)*x.^2+0.26340.*x+0.017839;fy=abs(f-y);fy2=fy.^2;Ew=max(fy),E1=sum(fy)/n,E2=sqrt((sum(fy2))/n)plot(xi,y,'t*'), hold on, plot(t,F,'b-'), hold off所得函数为4332(t)0.3436410 5.2156100.26340.013839t t t ϕ--=⨯-⨯++ 运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差Ew ，平均误差E1和均方根误差E2及其数据点),(i i y x 和拟合曲线y=f(x)的图形如图5.1.Ew =0.4243E1 =0.0911E2 =0.1467图5.1一元三次多项式拟合曲线误差图5.2用一元四次多项式()4433221t a t a t a t a t +++=ϕ进行拟合：计算多项式系数：a 1, a 2, a 3, a 4xi=[0 5 10 15 20 25 30 35 40 45 50 55];y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64];n=length(xi);x=0:0.01:55;f=0.6026.*10.^(-6)*x.^4-0.31918.*10.^(-4)*x.^3-0.0029323.*x.^2+0.23807.*x+0.060449;x=0:0.01:55;F=0.6026.*10.^(-6)*x.^4-0.31918.*10.^(-4)*x.^3-0.0029323.*x.^2+0.23807.*x+0.060449;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所得函数为644332(t)0.6026100.3191810 2.9323100.238070.060449t t t t ϕ---=⨯-⨯-⨯++运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差Ew ，平均误差E1和均方根误差E2及其数据点),(i i y x 和拟合曲线y=f(x)的图形如图5.2。
Ew = 0.3897 E1 = 0.1034、 E2 =0.1429图5.2一元四次多项式拟合曲线误差图5.3用一元二次多项式()2321t a t a a t ++=ϕ进行拟合：计算多项式系数：a 1, a 2, a 3输入程序： >> syms a1 a2 a3x=[0 5 10 15 20 25 30 35 40 45 50 55]; fi=a1.*x.^2+ a2.*x+ a3运行后屏幕显示关于a 1,a 2和a 3的线性方程组：fi=[ a3, 25*a1 + 5*a2 + a3, 100*a1 + 10*a2 + a3, 225*a1 + 15*a2 + a3, 400*a1 + 20*a2 + a3, 625*a1 + 25*a2 + a3, 900*a1 + 30*a2 + a3, 1225*a1 + 35*a2 + a3, 1600*a1 + 40*a2 + a3, 2025*a1 + 45*a2 + a3, 2500*a1 + 50*a2 + a3, 3025*a1 + 55*a2 + a3]编写构造误差平方和的MATLAB 程序：y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64];fi =[ a3, 25*a1 + 5*a2 + a3, 100*a1 + 10*a2 + a3, 225*a1 + 15*a2 + a3, 400*a1 + 20*a2 + a3, 625*a1 + 25*a2 + a3, 900*a1 + 30*a2 + a3, 1225*a1 + 35*a2 + a3, 1600*a1 + 40*a2 + a3, 2025*a1 + 45*a2 + a3, 2500*a1 + 50*a2 + a3, 3025*a1 + 55*a2 + a3];fy=fi-y;fy2=fy.^2;J=sum(fy.^2) 运行后屏幕显示误差平方和如下：J =(100*a1 + 10*a2 + a3 - 54/25)^2 + (25*a1 + 5*a2 + a3 - 127/100)^2 + (225*a1 + 15*a2 + a3 - 143/50)^2 + (400*a1 + 20*a2 + a3 - 86/25)^2 + (900*a1 + 30*a2 + a3 - 83/20)^2 + (625*a1 + 25*a2 + a3 - 387/100)^2 + (1225*a1 + 35*a2 + a3 - 437/100)^2 + (1600*a1 + 40*a2 + a3 - 451/100)^2 + (2025*a1 + 45*a2 + a3 - 229/50)^2 + (2500*a1 + 50*a2 + a3 - 201/50)^2 + (3025*a1 + 55*a2 + a3 - 116/25)^2 + a3^2为求4321,,,a a a a 使J 达到最小，只需利用极值的必要条件0=∂∂ka J)4,3,2,1(=k ，得到关于4321,,,a a a a 的线性方程组，这可以由下面的MATLAB 程序完成，即输入程序 ：>> syms a1 a2 a3J =(100*a1 + 10*a2 + a3 - 54/25)^2 + (25*a1 + 5*a2 + a3 - 127/100)^2 + (225*a1 + 15*a2 + a3 - 143/50)^2 + (400*a1 + 20*a2 + a3 - 86/25)^2 + (900*a1 + 30*a2 + a3 - 83/20)^2 + (625*a1 + 25*a2 + a3 - 387/100)^2 + (1225*a1 + 35*a2 + a3 - 437/100)^2 + (1600*a1 + 40*a2 + a3 - 451/100)^2 + (2025*a1 + 45*a2 + a3 - 229/50)^2 + (2500*a1 + 50*a2 + a3 - 201/50)^2 + (3025*a1 + 55*a2 + a3 - 116/25)^2 + a3^2;Ja1=diff(J,a1);Ja2=diff(J,a2);Ja3=diff(J,a3);Ja11=simple(Ja1),Ja21=simple(Ja2),Ja31=simple(Ja3), 运行后屏幕显示J 分别对a1, a2 ,a3的偏导数如下： Ja11 =49967500*a1 + 1089000*a2 + 25300*a3 - 217403/2 Ja21 = 1089000*a1 + 25300*a2 + 660*a3 - 27131/10 Ja31 = 25300*a1 + 660*a2 + 24*a3 - 3987/50解线性方程组Ja11 =0，Ja21 =0，Ja31 =0输入下列程序： >>A=[49967500,1089000,25300;1089000, 25300,660;25300,660,24]; B=[217403/2,27131/10,3987/50]; C=B/A, F=poly2sym(C)运行后屏幕显示拟合函数f 及其系数C 如下： C =-0.0024 0.2037 0.2305F = (7338734818964133*x)/36028797018963968 - (5489104202452799*x^2)/2305843009213693952 + 8303449950332545/36028797018963968故所求的拟合曲线为：0.230470.203690.0023805)(22++-=x x x f编写下面的MATLAB 程序估计其误差，并作出拟合曲线和数据的图形。