}//主函数结束
//*
******************************************************************
doublenewtonz（doublex）{
doublex1=0.0,t;
//牛顿迭代子函数
t=（7*pow（x,6）-4*28*pow（x,3））;
if（t==0）
exit（0）;
x1=（x-（（pow（x,7）-28*pow（x,4）+14）/t））;
returnx1;
}
double*newton（doublea,doubleb,doubleeps）{
doublex0=0.0,x1=1.0,x2=0.0,re[2];intk=0;
//牛顿迭代函数
x0=a;
while（x0>eps）
{
k++;
x2=x1;
x1=newtonz（x1）;x0=fabs（x1-x2）;
}re[0]=x1;
//代入a迭代计算
//调用牛顿迭代子函数
x0=b,k=0;
while（x0>eps）
{
k++;
x2=x1;
x1=newtonz（x1）;
x0=fabs（x1-x2）;
}re[1]=x1;
//代入b迭代计算
//调用牛顿迭代子函数
returnre;
1.3计算结果打印
■1'C:
\U5er5\AdminiStratoDesl'lop、.D亡bug\胡亚梅eye"
li.845497
Pi'essanj/ke^tocontinue
1.4MATLAB上机程序
functiony=Newton（f,df,xO,eps,M）d=0;
fork=1:
M
iffeval（df,xO）==O
d=2;break
else
x1=xO-feval（f,xO）/feval（df,xO）;
end
e=abs（x1-x0）;
x0=x1;
ife<=eps&&abs（feval（f,x1））v=epsd=1;break
end
end
ifd==1
y=x1;
elseifd==0
y='迭代M次失败：
else
y='奇异’
end
functiony=df（x）
y=7*xA6-28*4*xA3;
End
functiony=f（x）
y=xA7-28*xA4+14;
End
>>x0=1.9;
>>eps=0.00001;
>>M=100;
>>x=Newton（'f,'df,x0,eps,M）;
>>vpa（x,7）
1.5问题讨论
1.需注意的是，要使用Newton迭代法须f（x）x728x414满足定理中的条件I,n,川，以及
f（xc）•f''（x°）>0。
要用误差范围来控制循环的次数，保证循环的次数和质量，编写程序过程中要
注意标点符号的使用，正确运用适当的标点符号，Newton迭代法是局部收敛的，在使用时应先确定
初始值，否则所得的解可能不在所要求的范围内
（3）因为newton法求方程是平方收敛的，所以较为精确，但是要求出函数的导数，且必须有二阶导数。
第二题
2.已知函数值如下表：
X
1
2
3
4
5
f（x）
0
0.69314718:
1.0986123
1.3862944
1.6094378
X
6
7
8
9
10
f（x）
1.7917595
1.9459101
2.079445
2.1972246
2.3025851
f'（x）
f'
（1）=1
f'（10）=0.1
试用三次样条插值求f（4.563）及f'（4.563）的近似值
2.1理论依据及方法应用条件三次样条插值函数可定义为：
对于［a,b］上的一个划分n
a=2）
如果定义在［a,b］上的函数S（x）,满足
（1）.在［xi,xi+i］上为3次多项式；
（2）.S（x）,S'（x）,
S"（x）在［a,b］上连续，则称S（x）在［a,b］上划分的3次样条函数，如果对于f（x）［a,b］,
s（x）f（x）还满足s（xjf（x）,i0~n，则称s（x）为f（x）的三次样条插值函数。
其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中1,x,,x2,x3,（x-xj）+3为基函数,而取B样条函数Q3（（x-Xj）/h）为基函数.由于三次样条函数空间是N+3隹的,故我们把分点扩大到X1,Xn+1则任意三次样条函数可用Q3（（X-Xj）/h）线性组合来表示S（x）=N1CjQ3（（x-xj）/h）这样对不同插值
j1
问题,若能确定Cj由解的唯一性就能求得S（x）。
由s（xi）=yi,I=1,2,…Ns'（x0）=y0,s（xN）=yn可得
N1
S（Xi）=CjQ3（（Xi-Xj）/h）=yi
（x°）=1/h
CjQ3
（X0-Xj）
/h）=y
（xN）=1/h
CjQ3
（XN-Xj）
/h）=y
d0
d1
dN
d0
A—y'o）
00
d0
yy.
（j1jh•4h.（
j1j
h.
j
（yN
hN1
yN
yj
hj1
yj1
）6f（Xj1，Xj,Xj1）
/*宏定义*/
/*追赶法求解三弯矩方程*/
hj1
hj
2.2计算程序
#include
#include
#defineN10
main（）
{
floats,ds,t;
floatdy0=1,dy9=0.1;
intj;
intx[N]={1,2,3,4,5,6,7,8,9,10};
floaty[N]={0,0.69314718,1.0986123,1.3862944,1.6094378,
1.7917595,1.9459101,2.079445,2.1972246,2.3025851};
intb[N]={2,2,2,2,2,2,2,2,2,2},h[N-1];
floatd[N],u[N-1],v[N-1],a[N-1],c[N-1],B[N],l[N-1],p[N],X[N];
for（j=1;jv=9;j++）
h[j-1]=x[j]-x[j-1];
d[0]=6/h[0]*（y[1]/h[0]-y[0]/h[0]-dy0）;
d[9]=6/h[8]*（dy9-y[9]/h[8]+y[8]/h[8]）;
for（j=1;j<=8;j++）d[j]=6/（h[j-1]+h[j]）*（y[j+1]/h[j]-y[j]/h[j]-y[j]/h[j-1]+y[j-1]/h[j-1]）;
for（u[8]=1,j=0;j<=7;j++）
u[j]=h[j-1]/（h[j-1]+h[j]）;
for（v[0]=1,j=1;jv=8;j++）
v[j]=h[j]/（h[j-1]+h[j]）;
for（j=0;j<=8;j++）
a[j]=u[j];
for（j=0;j<=8;j++）
c[j]=v[j];
for（B[0]=b[0],j=1;jv=9;j++）
B[j]=b[j]-a[j]/B[j-1]*c[j-1];
for（j=1;j<=9;j++）
l[j]=a[j]/B[j-1];
for（j=1;j<=9;j++）
p[j]=d[j]-l[j]*p[j-1];
X[9]=p[9]/B[9];
for（j=8;j>=0;j--）
X[j]=p[j]/B[j]-c[j]*X[j+1]/B[j];
t=4.563;
/*解f（x）的值
/*解f'（x）
/*打印结果*/
s=X[3]*pow（（x[4]-t）,3）/6/h[3]+X[4]*pow（（t-x[3]）,3）/6/h[3]+（y[3]-X[3]*h[3]*h[3]/6）*（x[4]/h[3]-t/h[3]）+（y[4]-X[4]*h[3]*h[3]/6）*（t/h[3]-x[3]/h[3]）;
*/
ds=-X[3]*pow（（x[4]-t）,2）/2/h[3]+X[4]*pow（（t-x[3]）,2）/2/h[3]-（y[3]-X[3]*h[3]*h[3]/6）/h[3]+（y[4]-X[4]*h[3]*h[3]/6）/h[3];的值*/
printf（"s=%f\nds=%f\n",s,ds）;
}2.3计算结果打印
2.4MATLAB上机程序
functionQ=san（ssss,p）
Q=zeros（2,1）;
x=[1;2;3;4;5;6;7;8;9;10];
y=[0;0.69314718;1.0986123;1.3862944;1.6094378;1.7917595;1.9459101;2.079445;2.1972246;2.3025851];h=zeros（10,1）;
d=zeros（10,1）;
u=zeros（10,1）;
v=zeros（10,1）;
r=zeros（10,1）;
l=zeros（10,1）;
z=zeros（10,1）;
m=zeros（10,1）;
fort=1:
1:
9;
h（t）=x（t+1）-x（t）;
end
d
（1）=6/h
（1）*（（y
（2）-y
（1））/h
（1）-1）;
d（10）=6/h（9）*（0.1-（y（10）-y（9））/h（9））;
fort=1:
1:
8
u（t+1）=h（t）/（h（t）+h（t+1））;
v（t+1）=1-u（t+1）;
d（t+1）=6/（h（t）+h（t+1））*（（y（t+2）-y（t+1））/（x（t+2）-x（t+1））-（y（t+1）-y（t））/（x（t+1）-x（t）））;
end
u（10）=1;v
（1）=1;r
（1）=d
（1）;
fort=2:
1:
10
l（t）=u（t）/r（t-1）;
r（t）=d（t）-l（t）*v（t-1）;
end
z
（1）=d
（1）;
fort=2:
1:
10
z（t）=d（t）-l（t）*z（t-1）;
end
m（10）=z（10）/r（10）;
fort=9:
-1:
1
m（t）=（z（t）-v（t）*m（t+1））/r（t）;
end
fort=1:
1:
10
ifp>=t&&p<（t+1）
Q（:
1）=feval（ssss,p,t,x,m,h,y）;break
end
end
functionQ=ssss（p,t,x,m,h,y）
Q=zeros（2,1）;
Q（1,1）=（（power（（x（t+1）-p）,3）*m（t）+power（（p-x（t））,3）*m（t+1））/6+（y（t）-m（t）*h（t）*h（t）/6）*（x（t+1）-p）+（y（t+1）-m（t+1）*h（t）*h（t）/6）*（p-x（t）））/h（t）;
Q（2,1）=（-（power（（x（t+1）-p）,2）*m（t）+power（（p-x（t））,2）*m（t+1））/2+（y（t）-m（t）*h（t）*h（t）/6）+（y（t+1）-m（t+1）*h（t）*h（t）/6））/h（t）;end
2.5问题讨论
1.若要用追赶法求解三对角方程组，三对角阵需要满足：
A（i=1,2,…,n）均非奇异，保证A有唯一的
Doolittle分解；aG丰0;
2.样条插值效果比Lagrange插值好,三次样条插值的解存在且唯一，近似误差较小.并且没有Runge现象。
第二题
3.用Romberg法求：
3
xx1.4（5x7）sinx2dx（允许误差&
3.1理论依据及方法应用条件弹b?
f（a）f（b）数值积分的Romberg算法计算步骤如下：
⑴1（I1）
T1T1
2
4mTmk）Tm
=0.00001）。
（0）（0）
111
时,
（k1）
Tm1
就停机
ba
（k1）
211
fa
i1
（2i
m1,2,,lk1,2,,lm1
3.2计算程序
/*定义函数f（x）*/
#include#include#defineN9floatf（floatx）
{
floaty;
y=pow（3,x）*pow（x,1.4）*（5*x+7）*sin（x*x）;
return（y）;
}
main（）
{
floatT[N+1][N+1],h[N+1],a=1,b=3,m[N+1];
inti,l;
T[1][0]=（b-a）*（f（a）+f（b））/2;
l=1;
while（l<=N）
{
m[l]=0;
for（i=1;i<=（pow（2,l-1））;i++）
m[l]+=f（a+（2*i-1）*（b-a）/pow（2,l））;
T[1][l]=（T[1][l-1]+（b-a）*m[l]/pow（2,l-1））/2;
l++;
}
i=1;
while（i<=N）
{
for（l=1;l<=N-i+1;l++）
T[i+1][l-1]=（pow（4,i）*T[i][l]-T[i][l-1]）/（pow（4,i）-1）;
h[i]=T[i][0]-T[i+1][0];
if（fabs（h[i]）<=1e-5）break;
i++;
}
printf（"Theansweris:
%f",T[i+1][0]）;
}
3.3计算结果打印
3.4MATLAB上机程序
function[T,n]=mromb（f,a,b,eps）
ifnarginv4,eps=1e-6;end
h=b-a;
R（1,1）=（h/2）*（feval（f,a）+feval（f,b））;
n=1;J=0;err=1;
while（err>eps）
J=J+1;h=h/2;S=0;
fori=1:
n
x=a+h*（2*i-1）;
S=S+feval（f,x）;
end
R（J+1,1）=R（J,1）/2+h*S;
fork=1:
J
R（J+1,k+1）=（4Ak*R（J+1,k）-R（J,k））/（4Ak-1）;
end
err=abs（R（J+1,J+1）-R（J+1,J））;
n=2*n;
end
R;
T=R（J+1,J+1）;
formatlong
f=@（x）（3.Ax）*（x.A1.4）*（5*x+7）*sin（x*x）;
[T,n]=mromb（f,1,3,1.e-5）
3.5问题讨论
1、Romberge算法的优点是：
a、把积分化为代数运算，而实际上只需求「⑴,以后用递推可得
b、算法简单且收敛速度快，一般4或5次即能达到要求。
c、节省存储量，算出的可存入。
2、Romberge算法的缺点是：
a、对函数的光滑性要求较高。
b、计算新分点的值时，这些数值的个数成倍增加。
第四题
4.用定步长四阶Runge-Kutta法求解
厂dy,/dt1
dyz/dty3
dy3/dt10001000y2100y3
yi（0）0
y2（0）0
Iy3（0）0
h0.0005,打印yi（0.025）,y（0.045）,%（0.085）,%（0.1）,（i1,2,3）
4.1理论依据及方法应用条件
Runge-Kutta法的基本思想：
旳1不是按Taylor公式展开，而是先写成tn处附近的值的线性组合（有待定系数）,再按Taylor公式展开，然后确定待定常数。
四阶古典Runge-Kutta公式：
Yn1
Yn（K12K
6
22K3K4）
K1
hF（Xn,Yn）
K2
hF（Xn
£h,Yn
2
2K1）
2
K3
hF（Xn
^h,Yn
2
丄心）
2
K4
hF（Xn
£h,Yn
2
1
-K3）
2
4.2计算程序
#includeintmain（）
{
inti;
doubleh=0.0005;
doublek1,k2,k3,k4;
doubley1=0.0,y2=0.0,y3=0.0;
for（i=1;i<=200;i++）
{
k1=k2=k3=k4=h*1.0;
y1+=（k1+2*k2+2*k3+k4）/6;
k仁k2=k3=k4=h*y3;y2+=（k1+2*k2+2*k3+k4）/6;
k仁h*（1000-1000*y2-100*y3）;k2=h*（1000-1000*y2-100*（y3+0.5*k1））;k3=h*（1000-1000*y2-100*（y3+0.5*k2））;k4=h*（1000-1000*y2-100*（y3+k3））;y3+=（k1+2*k2+2*k3+k4）/6;
if（i==50）
{
printf（"\ny1（0.025）=%fcontinue;
}
if（i==90）
{
printf（"\ny1（0.045）=%fcontinue;
}
if（i==170）
{
printf（"\ny1（0.085）=%fcontinue;
}
if（i==200）
printf（"\ny1（0.100）=%f
}
y2（0.025）=%f
y2（0.045）=%f
y2（0.085）=%f
y2（0.100）=%f
y3（0.025）=%f",y1,y2,y3）;
y3（0.045）=%f",y1,y2,y3）;
y3（0.085）=%f",y1,y2,y3）;
y3（0.100）=%f\n\n",y1,y2,y3）;
4.3计算结果打印
4.4MATLAB上机程序
functionY=R_K（df1,a,b,h）
m=（b-a）/h;
Y=zeros（3,1）;
S=zeros（3,1）;
K=zeros（3,4）;
x=a;y1=a;y2=a;y3=a;
forn=1:
m
K（:
1）=feval（df1,x,y1,y2,y3）;
x=x+0.5*h;S（:
1）=Y+0.5*h.*K（:
1）;
y1=S（1,1）;y2=S（2,1）;y3=S（3,1）;
K（:
2）=feval（df1,x,y1,y2,y3）;
S（:
1）=Y+0.5*h.*K（:
2）;
y1=S（1,1）;y2=S（2,1）;y3=S（3,1）;
K（:
3）=feval（df1,x,y1,y2,y3）;
x=x+0.5*h;S（:
1）=Y+h.*K（:
3）;
y1=S（1,1）;y2=S（2,1）;y3=S（3,1）;
K（:
4）=feval（df1,x,y1,y2,y3）;
Y=Y+h.*（K（:
1）+2.*K（:
2）+2.*K（:
3）+K（:
4））/6;
end
functionZ=df1（x,y1,y2,y3）
Z=zeros（3,1）;
Z
（1）=0*x+0*y1+0*y2+0*y3+1;
Z
（2）=0*x+0*y1+0*y2+1*y3;
Z（3）=0*x+0*y1-1000*y2-100*y3+1000;
end
4.5问题讨论
1.定步长四阶runge-kutta法稳定，精度高，可根据有y'f（t,y）变化的情况与需要的精度自动修
改步长，误差小且程序简单，存储量少
2.但是Runge-Kutta法需要每步都计算函数值f（t,y）四次,在函数较复杂时，工作量就会变得较大
可靠性有待核查。
第五题
5.已知A与b
升级会员 