数值分析第五版计算实习题复习进程
数值分析计算实习题二
《数值分析》计算实习题二算法设计方案1．主要计算步骤：计算函数f(x,y)在拟合所需的节点处的函数值。
将各拟合节点（x i，y j）分别带入非线性方程组0.5 cos t + u + v + w – x = 2.67t + 0.5 sin u + v + w – y = 1.070.5t + u + cos v + w – x =3.74t + 0.5u + v + sin w – y =0.79解非线性方程组得解向量（t ij，u ij，v ij，w ij）。
对数表z(t,u)进行分片二次代数插值，求得对应（t ij，u ij）处的值，即为f(x i，y j) 的值。
对上述拟合节点分别进行x，y最高次数为k（k=0,1,2,3…）次的多项式拟合。
每次拟合后验证误差大小，直到满足要求。
2．求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。
3．对z(t，u)进行插值选择分片二次插值。
4．拟合基函数φr(x)ψs(y)选择为φr(x)=x r，ψs(y)=y s。
拟合系数矩阵c通过连续两次解线性方程组求得。
一．源程序#include "stdio.h"#include "stdlib.h"#include "math.h"void Doolittle(double *A,int n,int *M)//功能说明：对n阶矩阵A进行选主元的Doolittle分解//参数说明：A：欲进行分解的方阵，同时也是返回参数，分解后的结果// 存储于A中// n：方阵A的维数// M；（返回参数）n维向量，记录选主元过程中行交换的次序{int i,j,k,t;double *s;double Maxs,temp;s=(double*) calloc(n,sizeof(double));for(k=0;k<n;k++){for(i=k;i<n;i++){s[i]=A[i*n+k];for(t=0;t<k;t++) s[i]-= A[i*n+t] * A[t*n+k];}Maxs=abs(s[k]); M[k]=k;for(i=k+1;i<n;i++){if(Maxs<abs(s[i])){Maxs=abs(s[i]);M[k]=i;}}if(M[k]!=k){for(t=0;t<n;t++){temp=A[k*n+t];A[k*n+t]=A[M[k]*n+t];A[M[k]*n+t]=temp;}temp=s[k];s[k]=s[M[k]];s[M[k]]=temp;}if(Maxs<(1e-14)){s[k]=1e-14;printf("%.16e方阵奇异\n",Maxs);}A[k*n+k]=s[k];for(j=k+1;(j<n)&&(k<n-1);j++){for(t=0;t<k;t++) A[k*n+j]-=A[k*n+t]*A[t*n+j];A[j*n+k]=s[j]/A[k*n+k];}}}void Solve_LUEquation(double* A,int n,double* b,double* x) //功能说明：解方程LUx＝b，其中L、U共同存储在A中//参数说明：A：经Doolittle分解后的方阵// n：方阵A的维数// b：方程组的右端向量// x：（返回参数）方程组的解向量{int i,t;for(i=0;i<n;i++){x[i]=b[i];for(t=0;t<i;t++) x[i]-=A[i*n+t]*x[t];}for(i=n-1;i>-1;i--){for(t=i+1;t<n;t++) x[i]-=A[i*n+t]*x[t];x[i]/=A[i*n+i];}}void Transpose(double *A,int m,int n,double* AT)//功能说明：求m×n阶矩阵A的转置AT//参数说明：A：已知m×n阶矩阵// m：A的行数// n：A的列数// AT：（返回参数）A的转置矩阵（n×m）{int i,j;for(i=0;i<m;i++)for(j=0;j<n;j++) AT[j*m+i]=A[i*n+j];}void Solve_LEquation(double* A,int n,double* B,double* x,int m) //功能说明：解线性方程组Ax＝B，该函数可对系数矩阵相同// 而右端向量不同的多个方程组同时求解。
数值分析计算实习题
数值分析计算实习题-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN《数值分析》计算实习题姓名：学号：班级：第二章1、程序代码Clear;clc;x1=[ ];y1=[ ];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms x df d;df(1)=1;d(1)=y1(1);for i=2:n %求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i-1)*df(i);endP4=vpa(sum(d),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, 'variational');%调用三次样条函数q=;q1=q(1,:)*[^3;^2;;1];q1=vpa(collect(q1),5)q2=q(1,:)*[^3;^2;;1];q2=vpa(collect(q2),5)q3=q(1,:)*[^3;^2;;1];q3=vpa(collect(q3),5)q4=q(1,:)*[^3;^2;;1];q4=vpa(collect(q4),5)%求解并化简多项式2、运行结果P4 =*x - *(x - *(x - - *(x - *(x - *(x - - *(x - *(x - *(x - *(x - + q1 =- *x^3 + *x^2 - *x +q2 =- *x^3 + *x^2 - *x + q3 =- *x^3 + *x^2 - *x + q4 =- *x^3 + *x^2 - *x +3、问题结果4次牛顿差值多项式4()P x = *x - *(x - *(x - - *(x - *(x - *(x - - *(x - *(x - *(x - *(x - +三次样条差值多项式()Q x0.10.20.30.40.50.60.70.80.910.40.50.60.70.80.911.1323232321.33930.803570.40714 1.04,[0.2,0.4]1.3393 1.60710.88929 1.1643,[0.4,0.6]1.3393 2.4107 1.6929 1.4171,[0.6,0.8]1.3393 3.21432.8179 1.8629,[0.8,1.0]x x x x x x x x x x x x x x x x ⎧-+-+∈⎪-+-+∈⎪⎨-+-+∈⎪⎪-+-+∈⎩第三章1、程序代码Clear;clc; x=[0 1]; y=[1 ];p1=polyfit(x,y,3)%三次多项式拟合 p2=polyfit(x,y,4)%四次多项式拟合 y1=polyval(p1,x);y2=polyval(p2,x);%多项式求值plot(x,y,'c--',x,y1,'r:',x,y2,'y-.')p3=polyfit(x,y,2)%观察图像，类似抛物线，故用二次多项式拟合。
数值分析第五版第5章学习资料
n
即 de(A t) aijAij (i1,2,,n), j1
其中 A ij 为 a ij 的代数余子式，Aij(1)ijMij, M ij 为元素 a ij 的余子式.
行列式性质：
( ad ) ( A e ) d t B ( A e )d ( t B )A e , ,B t R n n .
有非零解，故系数行列式 deIt (A)0，记
a11 a12 p()det(I A) a21 a22
a1n a2n
（1.3）
an1 an2 ann n c1n1cn1cn 0.
p()称为矩阵 A的特征多项式，方程(1.3)称为矩阵 A的特
征方程.
9
因为 n次代数方程 p() 在复数域中有 n个根
其中用 ri 表示矩阵的第 i行. 由此看出，用消去法解方程组的基本思想是用逐次消
去未知数的方法把原方程组 Axb化为与其等价的三角 形方程组，而求解三角形方程组可用回代的方法.
上述过程就是用行的初等变换将原方程组系数矩阵化 为简单形式(上三角矩阵)，从而将求解原方程组(2.1)的 问题转化为求解简单方程组的问题.
n
n
trA aii i.
i1
i1
（1.4） （1.5）
称 trA为 A的迹.
A的特征值 和特征向量 x还有一下性质：
(1) AT 与 A有相同的特征值 及特征向量 .
(2)若 A非奇异，则 A1 的特征值为 1，特征向量为 x.
(3)相似矩阵 BS1AS有相同的特征多项式.
11
例1 求 A的特征值及谱半径
4x2x3 5,
2x3 6.
显然，方程组(2.6)是容易求解的，解为
x (1,2,3)T.
数值分析(第五版)计算实习题第五章作业教学资料
>> format compact
>> A=[3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34];
>> b=[1;1;1];
>> [RA,RB,n,X]=liezhu(A,b),h=det(A),C=cond(A)
输出：
请注意：因为RA=RB，所以方程组有唯一解
ans =
-9.5863 18.3741 -3.2258 3.5240
xX =
10.4661
jxX =
0.9842
Xgxx =
22.7396
xAb =
0.0076
xAbj =
0.0076
Acp =
2.9841e+03
第四题：
（1）输入：
建立m文件:
forn=2:6
a=hilb(n);
pnH(n-1)=cond(a,inf);
RA =
3
RB =
3
n =
3
X =
1.0e+03 *
1.5926
-0.6319
-0.4936
h =
-0.0305
C =
3.0697e+04
(2)输入：
>> A=[3.00 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34];
>> b=[1;1;1];
>> [RA,RB,n,X]=liezhu(A,b),h=det(A)
>> r=b-H*X,deltax=X-x
输出：
X =
数值分析-复习及习题选讲
5、线性方程组的数值解法
1.了解Gauss消元法的基本思想,知道适用范围 顺序Gauss消元法:矩阵A的各阶顺序主子式都不为零. 主元Gauss消元法:矩阵A的行列式不为零. 2.掌握矩阵的直接三角分解法。
会对矩阵进行Doolittle分解(LU)、Crout分解及Cholesky分解。
熟练掌握用三角分解法求方程组的解。 了解平方根法和追赶法的思想。 3.了解向量和矩阵的范数的定义,会判定范数(三要素非负性、齐 次性、三角不等式);会计算几个常用的向量和矩阵的范数； 了解范数的等价性和向量矩阵极限的概念。 4.了解方程组的性态，会计算简单矩阵的条件数。
k n
f
( n 1)
(2)记(t)=(t-x)k,则yj=(xj)=(xj-x)k, j=0,1,…,n.于是
n ( t ) k (t x) k f (t ) y j l j (t ) n 1 (t ) ( x j x) l j (t ) j 0 j 0 (n 1)! 取t=x,则有 n ( x j x) k l j ( x) 0
收敛于(x)在I上的唯一不动点x*.
都收敛于方程的唯一根x*.
推论 若(x)在x*附近具有一阶连续导数,且|(x*)|<1, 则对充分接近 x*的初值x0,迭代法xk+1=(xk)收敛. 3. 了解迭代法收敛阶的概念,会求迭代法收敛的阶.了解Aitken加速 技巧.
xk 1 C (1) xkp阶收敛于x*是指: lim k x p k
7.设(x)=x4+2x3+5, 在区间[-3,2]上, 对节点x0= -3, x1=-1,求出(x)的
三次Hermite插值多项式在区间[x0,x1]上的表达式及误差公式.
数值分析第五版计算实习题
弟二草插值法3.卜列数据点的插值可以得到平方根函数的近似，在区间064］上作图。
(1〉用这9个点做8次多项式插值Q x)。
(2)用三次样条(第一边界条件)程岸求S(X)。
从得到结果石在［0.64］ 1：・哪个插值更粘确：在区间［0,1］ I：•两种插值哪个更精确?(1) 8次多项式插值：(1)8次多项式插值:首先建立新的M-file:输入如卜代码(此为拉格朗口插值的功能函数)并保存function f=Language(x,y,x0)%求Li知数据点的拉格朗Fl插值多项式%己知数据点的x坐标向量：x%已知数据点的y坐标向量:y%插值的x坐标：x0%求得的拉格朗H插值多项式或在X0处的插值：fsyms t;ifi(lcngth(x)=length(y))n=length(x);elsedisp(*x和y的维数不相等!)；return;end %检错tbr(i=l:n)i=y(i);fbr(j=1:i-l)l=l*(t-x(j))/(x(i)-x(j))；end;for(j=i-M:n)end;for(j=i+l:n) l=l*(t-x(j))/(x(i)-x(j))； end;simplify(f);if(i==n) if|nargin=3)f=subs(C't\xO);else f=collcct(f);f=vpa(f,6);endendend再建立新的M-file：输入：clear;x=[0 1 49 16 25 36 49 64];y=[0：l：8]；%计算拉格朗口基丞数%计算拉格朗ri插值函数%化简%计算插值点的曲数值%将插值多项式展开%将插值多项式的系数化成6位精度的小数f=Uinguage(x,y) 运行得到f=1.32574*1-381410*t A2+.604294e-1 *t A3+.222972e-3 *t A5-.542921 e-5*t A6+.671268e・7T7・.328063e・9T8・.498071 e-2*t A4 这就是8次多项式插值L s(x)= 1.32574怜.381410*t A2+.604294e-1 *t A3+.222972e-3 *t A5-.542921 e-5*t A6+.671268e-7*t A7-.328063e-9*t A8-.498071 e-2*t A4. (2)三次样条插值：建立新的M-filc：输入：clear;x=[0 I 49 1625 36 4964];尸[0:8];t=[0:0.1:64];Y=t.A(0.5);O=Language(x,y)f= 1,32574*t-.381410*t.A2+.604294e-1 *t.A3+.222972e-3*t.A5-.542921 e・5*(. W+.671268e-7*t.A7-.328063e-9*t.A8-.498071 e-2 *t.A4;S=interp l(x,y,t.'spline,);plol(x,y,o；(・YY.lf.'b'」S'g:');grid;运行程序得到如下图：从结果屮很明显可以看出在[0.64].上.三次样条插值更精确，儿乎与原函数帀合。
数值分析第五版计算实习题
数值分析计算实习题第二章2-1程序：clear;clc;x1=[0.2 0.4 0.6 0.8 1.0];y1=[0.98 0.92 0.81 0.64 0.38];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms x df d;df(1)=1;d(1)=y1(1);for i=2:n %求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i)*df(i);enddisp('4次牛顿插值多项式');P4=vpa(collect((sum(d))),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, 'variational');%调用三次样条函数q=pp.coefs;disp('三次样条函数');for i=1:4S=q(i,:)*[(x-x1(i))^3;(x-x1(i))^2;(x-x1(i));1];S=vpa(collect(S),5)endx2=0.2:0.08:1.08;dot=[1 2 11 12];figureezplot(P4,[0.2,1.08]);hold ony2=fnval(pp,x2);x=x2(dot);y3=eval(P4);y4=fnval(pp,x2(dot));plot(x2,y2,'r',x2(dot),y3,'b*',x2(dot),y4,'co');title('4次牛顿插值及三次样条');结果如下：4次牛顿插值多项式P4 = - 0.52083*x^4 + 0.83333*x^3 - 1.1042*x^2 + 0.19167*x + 0.98三次样条函数x∈[0.2,0.4]时， S = - 1.3393*x^3 + 0.80357*x^2 - 0.40714*x + 1.04 x∈[0.4,0.6]时，S = 0.44643*x^3 - 1.3393*x^2 + 0.45*x + 0.92571x∈[0.6,0.8]时，S = - 1.6964*x^3 + 2.5179*x^2 - 1.8643*x + 1.3886 x∈[0.8,1.0]时，S =2.5893*x^3 - 7.7679*x^2 + 6.3643*x - 0.80571输出图如下2-3（1）程序：clear;clc;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];%插值点n=length(y1);a=ones(n,2);a(:,2)=-x1';c=1;for i=1:nc=conv(c,a(i,:));endq=zeros(n,n);r=zeros(n,n+1);for i=1:n[q(i,:),r(i,:)]=deconv(c,a(i,:));%wn+1/(x-xk) end 三次样条插值曲线4次牛顿插值曲线Dw=zeros(1,n);for i=1:nDw(i)=y1(i)/polyval(q(i,:),x1(i));%系数endp=Dw*q;syms x L8;for i=1:nL8(i)=p(n-i+1)*x^(i-1);enddisp('8次拉格朗日插值');L8=vpa(collect((sum(L8))),5)xi=0:64;yi=polyval(p,xi);figureplot(xi,yi,x1,y1,'r*');hold ontitle('8次拉格朗日插值');结果如下：8次拉格朗日插值L8 =- 3.2806e-10*x^8 + 6.7127e-8*x^7 - 5.4292e-6*x^6 + 0.00022297*x^5 - 0.0049807*x^4 + 0.060429*x^3 - 0.38141*x^2 + 1.3257*x输出图如下：第五章4-1（3）程序：clc;clear;y= (x)sqrt(x).*log(x);a=0;b=1;tol=1e-4;p=quad(y,a,b,tol);fprintf('采用自适应辛普森积分结果为： %d \n', p); 结果如下：采用自适应辛普森积分结果为： -4.439756e-01第九章9-1(a)程序：clc;clear;a=1;b=2;%定义域h=0.05;%步长n=(b-a)/h;y0=1;%初值f= (x,y) 1/x^2-y/x;%微分函数Xn=linspace(a,b,n+1);%将定义域分为n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值%以下根据改进欧拉公式求解for i=1:nxn=Xn(i);xnn=Xn(i+1);yn=Yn(i);yp=yn+h*f(xn,yn);yc=yn+h*f(xnn,yp);yn=(yp+yc)/2;Yn(i+1)=yn;endXn=Yn;%以下根据经典四阶R-K法公式求解for i=1:nxn=Xn(i);yn=Yn(i);k1=f(xn,yn);k2=f(xn+h/2,yn+h/2*k1);k3=f(xn+h/2,yn+h/2*k2);k4=f(xn+h,yn+h*k3);yn=yn+h/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;enddisp(' 改进欧拉法四阶经典R-K法'); disp([Xn' Yn'])结果如下：改进欧拉法四阶经典R-K法1 10.99887 0.998850.99577 0.99780.99114 0.996940.98532 0.996340.97857 0.996030.97111 0.996060.96311 0.996450.9547 0.997230.94598 0.998410.93705 10.92798 1.0020.91883 1.00440.90964 1.00730.90045 1.01060.89129 1.01430.88218 1.01840.87315 1.02290.86421 1.02780.85538 1.03310.84665 1.0388（b）程序：clc;clear;a=0;b=1;%定义域H=[0.1 0.025 0.01];%步长y0=1/3;%初值f= (x,y) -50*y+50*x^2+2*x;%微分函数xi=linspace(a,b,11);Y=1/3*exp(-50*xi)+xi.^2;%准确解Ym=zeros(1,11);for j=1:3h=H(j);n=(b-a)/h;Xn=linspace(a,b,n+1);%将定义域分为n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值for i=1:nxn=Xn(i);yn=Yn(i);k1=f(xn,yn);k2=f(xn+h/2,yn+h/2*k1);k3=f(xn+h/2,yn+h/2*k2);k4=f(xn+h,yn+h*k3);yn=yn+h/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;endfor k=1:11m=0.1/h;Ym(k)=Yn(1+(k-1)*m);enddelta=Ym-Y;fprintf('步长为： %d \n', h);disp(' 四阶经典R-K法准确解误差'); disp([Ym' Y' delta'])end结果如下：步长为： 1.000000e-01四阶经典R-K法准确解误差0.33333 0.33333 04.6055 0.012246 4.593263.062 0.040015 63.022864.05 0.09 863.9611844 0.16 118431.6235e+05 0.25 1.6235e+052.2256e+06 0.36 2.2256e+063.0509e+07 0.49 3.0509e+074.1823e+08 0.64 4.1823e+085.7333e+09 0.81 5.7333e+097.8594e+10 1 7.8594e+10步长为： 2.500000e-02四阶经典R-K法准确解误差0.33333 0.33333 00.013015 0.012246 0.000768940.040063 0.040015 4.82e-050.090037 0.09 3.6857e-050.16004 0.16 3.6723e-050.25004 0.25 3.6722e-050.36004 0.36 3.6722e-050.49004 0.49 3.6722e-050.64004 0.64 3.6722e-050.81004 0.81 3.6722e-051 1 3.6722e-05步长为： 1.000000e-02四阶经典R-K法准确解误差0.33333 0.33333 00.012256 0.012246 9.5673e-060.040016 0.040015 7.8252e-070.090001 0.09 6.6347e-070.16 0.16 6.6226e-070.25 0.25 6.6225e-070.36 0.36 6.6225e-070.49 0.49 6.6225e-070.64 0.64 6.6225e-070.81 0.81 6.6225e-071 1 6.6225e-07由结果可知，步长越小，结果越精确。
李庆扬-数值分析第五版第5章和第7章习题答案解析复习课程
李庆扬-数值分析第五版第 5 章和第7 章习题答案解析第5章复习与思考题习题Ta iicii 1、设A 是对称阵且a n 0,经过高斯消去法一步后，A 约化为，证明A 2是对A 2称矩阵。
证明：所以A2为对称矩阵。
A (a j ))n i ；证明：(I )A 的对角元素—a H °(i I2|||,n) ； (2)A 2是对称正定矩阵;即A 是对称矩阵。
3、设L k 为指标为k 的初等下三角矩阵（除第~~k 列对角元以下元素外, L k 和单位阵I 相 同），即a ii °I2. 1.. a in设对称矩阵ai2A i2a22. …％ 2aina2 n. 1.. a nn a iiai2.. .ai2a22 —a i 2 ..an2A ⑴ aiiaiia22annaiia inan2a ^aina ii所以a Ta n2 —a i2..a ii[ai2...a n2]a nn — amaiiai23|ian2耳2 aiiai nannai n aiiOl n2、设A 是对称正定矩阵，经过高斯消去法一步后,A 约化为A (a j ) n ，其中A (O j )n ，(i )依次取 X i (0,0, ,0,i,0,所以有a iiX T A X0。
,0)T , i i,2,，n ，则因为A 是对称正定矩阵, (2) A 2中的元素满足a (2)aja i冋(i,jaii2,3, ,n),又因为A 是对称正定矩阵，满足a j a ji , i, ji,2,,n ，所以 a (2)aja ii ai jaiiajia iiaji(2)----- aji ■)a，则经过1次高斯校区法后，有0 a 2n亚 ％a l1a i nai111m k i,k 1m^k 1求证当i,j k时，L k I ij L k I ij也是一个指标为k的初等下三角矩阵，其中I j为初等置换矩阵。
4、试推导矩阵A的Crout分解A=LU的计算公式，其中L为下三角矩阵，U为单位上三角矩阵。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析计算实习题第二章2-1程序：clear;clc;x1=[0.2 0.4 0.6 0.8 1.0];y1=[0.98 0.92 0.81 0.64 0.38];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms x df d;df(1)=1;d(1)=y1(1);for i=2:n %求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i)*df(i);enddisp('4次牛顿插值多项式');P4=vpa(collect((sum(d))),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, 'variational');%调用三次样条函数q=pp.coefs;disp('三次样条函数');for i=1:4S=q(i,:)*[(x-x1(i))^3;(x-x1(i))^2;(x-x1(i));1];S=vpa(collect(S),5)endx2=0.2:0.08:1.08;dot=[1 2 11 12];figureezplot(P4,[0.2,1.08]);hold ony2=fnval(pp,x2);x=x2(dot);y3=eval(P4);y4=fnval(pp,x2(dot));plot(x2,y2,'r',x2(dot),y3,'b*',x2(dot),y4,'co');title('4次牛顿插值及三次样条');结果如下：4次牛顿插值多项式P4 = - 0.52083*x^4 + 0.83333*x^3 - 1.1042*x^2 + 0.19167*x + 0.98三次样条函数x∈[0.2,0.4]时， S = - 1.3393*x^3 + 0.80357*x^2 - 0.40714*x + 1.04 x∈[0.4,0.6]时，S = 0.44643*x^3 - 1.3393*x^2 + 0.45*x + 0.92571x∈[0.6,0.8]时，S = - 1.6964*x^3 + 2.5179*x^2 - 1.8643*x + 1.3886 x∈[0.8,1.0]时，S =2.5893*x^3 - 7.7679*x^2 + 6.3643*x - 0.80571输出图如下2-3（1）程序：clear;clc;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];%插值点n=length(y1);a=ones(n,2);a(:,2)=-x1';c=1;for i=1:nc=conv(c,a(i,:));endq=zeros(n,n);r=zeros(n,n+1);for i=1:n[q(i,:),r(i,:)]=deconv(c,a(i,:));%wn+1/(x-xk) end 三次样条插值曲线4次牛顿插值曲线Dw=zeros(1,n);for i=1:nDw(i)=y1(i)/polyval(q(i,:),x1(i));%系数endp=Dw*q;syms x L8;for i=1:nL8(i)=p(n-i+1)*x^(i-1);enddisp('8次拉格朗日插值');L8=vpa(collect((sum(L8))),5)xi=0:64;yi=polyval(p,xi);figureplot(xi,yi,x1,y1,'r*');hold ontitle('8次拉格朗日插值');结果如下：8次拉格朗日插值L8 =- 3.2806e-10*x^8 + 6.7127e-8*x^7 - 5.4292e-6*x^6 + 0.00022297*x^5 - 0.0049807*x^4 + 0.060429*x^3 - 0.38141*x^2 + 1.3257*x输出图如下：第五章4-1（3）程序：clc;clear;y= @(x) sqrt(x).*log(x);a=0;b=1;tol=1e-4;p=quad(y,a,b,tol);fprintf('采用自适应辛普森积分结果为： %d \n', p); 结果如下：采用自适应辛普森积分结果为： -4.439756e-01第九章9-1(a)程序：clc;clear;a=1;b=2;%定义域h=0.05;%步长n=(b-a)/h;y0=1;%初值f= @(x,y) 1/x^2-y/x;%微分函数Xn=linspace(a,b,n+1);%将定义域分为n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值%以下根据改进欧拉公式求解for i=1:nxn=Xn(i);xnn=Xn(i+1);yn=Yn(i);yp=yn+h*f(xn,yn);yc=yn+h*f(xnn,yp);yn=(yp+yc)/2;Yn(i+1)=yn;endXn=Yn;%以下根据经典四阶R-K法公式求解for i=1:nxn=Xn(i);yn=Yn(i);k1=f(xn,yn);k2=f(xn+h/2,yn+h/2*k1);k3=f(xn+h/2,yn+h/2*k2);k4=f(xn+h,yn+h*k3);yn=yn+h/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;enddisp(' 改进欧拉法四阶经典R-K法'); disp([Xn' Yn'])结果如下：改进欧拉法四阶经典R-K法1 10.99887 0.998850.99577 0.99780.99114 0.996940.98532 0.996340.97857 0.996030.97111 0.996060.96311 0.996450.9547 0.997230.94598 0.998410.93705 10.92798 1.0020.91883 1.00440.90964 1.00730.90045 1.01060.89129 1.01430.88218 1.01840.87315 1.02290.86421 1.02780.85538 1.03310.84665 1.0388（b）程序：clc;clear;a=0;b=1;%定义域H=[0.1 0.025 0.01];%步长y0=1/3;%初值f= @(x,y) -50*y+50*x^2+2*x;%微分函数xi=linspace(a,b,11);Y=1/3*exp(-50*xi)+xi.^2;%准确解Ym=zeros(1,11);for j=1:3h=H(j);n=(b-a)/h;Xn=linspace(a,b,n+1);%将定义域分为n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值for i=1:nxn=Xn(i);yn=Yn(i);k1=f(xn,yn);k2=f(xn+h/2,yn+h/2*k1);k3=f(xn+h/2,yn+h/2*k2);k4=f(xn+h,yn+h*k3);yn=yn+h/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;endfor k=1:11m=0.1/h;Ym(k)=Yn(1+(k-1)*m);enddelta=Ym-Y;fprintf('步长为： %d \n', h);disp(' 四阶经典R-K法准确解误差'); disp([Ym' Y' delta'])end结果如下：步长为： 1.000000e-01四阶经典R-K法准确解误差0.33333 0.33333 04.6055 0.012246 4.593263.062 0.040015 63.022864.05 0.09 863.9611844 0.16 118431.6235e+05 0.25 1.6235e+052.2256e+06 0.36 2.2256e+063.0509e+07 0.49 3.0509e+074.1823e+08 0.64 4.1823e+085.7333e+09 0.81 5.7333e+097.8594e+10 1 7.8594e+10步长为： 2.500000e-02四阶经典R-K法准确解误差0.33333 0.33333 00.013015 0.012246 0.000768940.040063 0.040015 4.82e-050.090037 0.09 3.6857e-050.16004 0.16 3.6723e-050.25004 0.25 3.6722e-050.36004 0.36 3.6722e-050.49004 0.49 3.6722e-050.64004 0.64 3.6722e-050.81004 0.81 3.6722e-051 1 3.6722e-05步长为： 1.000000e-02四阶经典R-K法准确解误差0.33333 0.33333 00.012256 0.012246 9.5673e-060.040016 0.040015 7.8252e-070.090001 0.09 6.6347e-070.16 0.16 6.6226e-070.25 0.25 6.6225e-070.36 0.36 6.6225e-070.49 0.49 6.6225e-070.64 0.64 6.6225e-070.81 0.81 6.6225e-071 1 6.6225e-07由结果可知，步长越小，结果越精确。
9-3主程序：clc;clear;options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-4]); [X,Y] = ode45(@rigid,[0 0.004],[1 1 0],options);subplot(1,2,1);plot(X,Y(:,1),'-',X,Y(:,2),'-.',X,Y(:,3),'.');title('四阶R-K方法');legend('y1', 'y2','y3');[X,YY] = ode23t(@rigid,[0 0.004],[1 1 0],options); subplot(1,2,2);plot(X,YY(:,1),'-',X,YY(:,2),'-.',X,YY(:,3),'.'); title('梯形法');legend('y1', 'y2','y3');rigid函数function dy = rigid (~,y)dy = zeros(3,1); % a column vectordy(1) = -0.013*y(1)-1000*y(1)*y(2);dy(2) = -2500*y(2) * y(3);dy(3) =-0.013*y(1)-1000*y(1)*y(2)-2500*y(2) * y(3);结果如图：。