大连理工大学矩阵与数值分析上机作业代码

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

《矩阵与数值分析》上机作业实验报告

陈晓学号：21204047机械工程学院教学班号：02任课教师：金光日

习题1

给定n阶方程组Ax=b，其中

61 7 861 15 A= ，b= 861 15 14 86

则方程组有解x=(1,1, ,1)。对n=10和n=84，分别用Gauss消去法和列主元消去法解方程组，并比较计算结果。T

1.1程序

（1）高斯消元法

n=10时，

>>[A,b]=CreateMatrix(10)

A=

6

8

0168000000001680000000016800000000168000000001680000000016800000000168000000001680000000016

b=

7

15

15

15

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

15

15

15

15

15

14

>>x=GaussMethod(A,b)

x=

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

n=84时，

>>[A,b]=CreateMatrix(84)

>>x=GaussMethod(A,b)

1.0e+008*

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

0.0000

……

-0.0001

0.0002

-0.0003

0.0007

-0.0013

0.0026

-0.0052

0.0105

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

-0.0209

0.0419

-0.0836

0.1665

-0.3303

0.6501

-1.2582

2.3487

-4.0263

5.3684

（2）列主元消去法

n=10时，

>>[A,b]=CreateMatrix(10)

>>x=MainElementMethod(A,b)

x=

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

n=84时，

>>[A,b]=CreateMatrix(84)

>>x=MainElementMethod(A,b)

x=

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

……

1.0000

1.0000

1.0000

1.0000

1.2计算结果分析

对于高斯消元法，当n=10时，得到的计算结果为x=（1.0000，1.0000，1.0000，1.0000，……，T1.0000），恰好为精确值；但当n=84时，得到的计算结果却为x=(0.0000,0.0000,0.0000,……，-0.0001,0.0002，……，2.3487，-4.0263，5.3684）T,完全偏离了精确值。

对于列主元消去法，当n=10时，得到的计算结果为x=（1.0000，1.0000，1.0000，

T1.0000，……，1.0000），恰好为精确值；当n=10时，得到的计算结果为x=（1.0000，1.0000，

1.0000，1.0000，……，1.0000）T，恰好为精确值。

由于使用高斯消元法时可能会出现小主元作除数的情况，虽然在求解系数矩阵为10阶的方程时，由于步数较少，尚能求得精确解，但在求解求解系数矩阵为84阶的巨型方程时，由于求解步数太多，使得小主元作除数的情况累加，最终导致“大数除小数”的情况发生，因而所得解严重偏离精确解。

而使用列主元消去法时，由于有效地避免了这一不利情况，因而无论是在求解系数矩阵为10阶还是84阶的方程时，都能得到精确解。

1.3M文件

.m1.3.1CreateMatrixCreateMatrix.m

function[A,b]=CreateMatrix(n)

%用于存放习题1的题目信息，并建立构造题目中矩阵的函数

%对矩阵A赋值

A1=6*ones(1,n);

A2=ones(1,n-1);

A3=8*ones(1,n-1);

A=diag(A1)+diag(A2,1)+diag(A3,-1);

%对向量b赋值

b=15*ones(n,1);

b(1)=7;

b(n)=14;

1.3.2GaussMethod.m

functionx=GaussMethod(A,b)

%高斯消元法返回方程Ax=b的解

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

%x是方程的解

%A为方程的系数矩阵

%b为方程的常数矩阵

n=length(b);

%定义单位下三角矩阵元胞数组，并将全部数组元素赋值为n阶单位阵

l=cell(1,n-1);

l(:)={eye(n)};

U=A;

fori=1:n-1

forj=i+1:n

l{i}(j,i)=-U(j,i)/U(i,i);

end

U=l{i}*U;

end

%计算L和U

L=inv(l{1});

fori=2:n-1

L=L/l{i};

end

x=U\(L\b);

.m1.3.3MainElementMethodMainElementMethod.m

functionx=MainElementMethod(A,b)

%列主元消去法返回方程Ax=b的解

%x是方程的解

%A为方程的系数矩阵

%b为方程的常数矩阵

n=length(b);

%定义单位下三角矩阵元胞数组，并将全部数组元素赋值为n阶单位阵

l=cell(1,n-1);

l(:)={eye(n)};

p=cell(1,n-1);

p(:)={eye(n)};

U=A;

fori=1:n-1

%找出列主元所在行，并交换首行，得到Pi

[Max,row]=max(U(i:n,i));

row=row+i-1;

mid=p{i}(row,:);

p{i}(row,:)=p{i}(i,:);

p{i}(i,:)=mid;

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

U=p{i}*U;

forj=i+1:n

l{i}(j,i)=-U(j,i)/Max;

end

U=l{i}*U;

end

%计算L

L=eye(n);

fori=1:n-2

forj=i+1:n-1

l{i}=p{j}*l{i}/p{j};

end

L=l{i}*L;

end

L=inv(L*l{n-1});

%计算P

P=eye(n);

fori=1:n-1

P=p{i}*P;

end

%计算x

x=U\(L\P*b);

习题2

设有方程组Ax=b，其中A∈R20×20，

0.5 0.25 3 3 0.5 0.5

0.25 0.5 A=

(a)选取不同的初始向量x(0) 0.5 0.25 0.53 0.5 0.25 0.53 和不同的右端向量b，给定迭代误差要求，用Jacobi和Gauss-Seidel迭代法计算，观测得出的迭代向量序列是否收敛。若收敛，记录迭代次数，分析计算结果并得出你的结论。

(b)选定初始向量初始向量x(0)和不同的右端向量b，如取x(0)=0,b=Ae,e=(1,1, 1)T。

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

将A的主对角线元素成倍增长若干次，非主对角元素不变，每次用Jacobi法计算，要求迭代误差满足x(k+1) x(k)

∞<10 6，比较收敛速度，分析现象并得出你的结论。

2.1程序

（1）取x(0)=，b=（1,1，……，1,1）T,迭代误差取默认值10 5，迭代次数上限取默认值50，使用Jacobi法进行迭代。

>>test2

>>b=ones(20,1)

>>x0=zeros(20,1)

>>[x,n]=JacobiMethod(A,b,x0)

x=

0.2766

0.2327

0.2159

0.2223

0.2227

0.2221

0.2222

0.2222

0.2222

0.2222

0.2222

0.2222

0.2222

0.2222

0.2221

0.2227

0.2223

0.2159

0.2327

0.2766

n=

17

（2）取x(0)=0，b=（1,1，……，1,1）T,迭代误差取默认值10 5，迭代次数上限取默认值50，使用Gauss-Seidel法进行迭代。

>>test2

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

>>b=ones(20,1)

>>x0=zeros(20,1)

>>[x,n]=Gauss_SeidelMethod(A,b,x0)

x=

0.2766

0.2327

0.2159

0.2223

0.2228

0.2221

0.2222

0.2222

0.2222

0.2222

0.2222

0.2222

0.2222

0.2222

0.2221

0.2228

0.2223

0.2159

0.2327

0.2766

n=

9

（3）取x(0)=0，b=（1,0,1,0,……，1,0）T,迭代误差取默认值10 5，迭代次数上限取默认值50，使用Jacobi法进行迭代。

>>test2

>>b=[1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0]'

>>x0=zeros(20,1)

>>[x,n]=JacobiMethod(A,b,x0)

x=

0.3238

-0.0985

0.3116

-0.0881

0.3110

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

0.3111

-0.0889

0.3111

-0.0889

0.3111

-0.0889

0.3111

-0.0889

0.3111

-0.0882

0.3105

-0.0957

0.3313

-0.0472

n=

16

（4）取x(0)=0，b=（1,0,1,0,……，1,0）T,迭代误差取默认值10 5，迭代次数上限取默认值50，使用Gauss-Seidel法进行迭代。

>>test2

>>b=[1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0]'

>>x0=zeros(20,1)

>>[x,n]=Gauss_SeidelMethod(A,b,x0)

x=

0.3238

-0.0985

0.3115

-0.0881

0.3110

-0.0889

0.3111

-0.0889

0.3111

-0.0889

0.3111

-0.0889

0.3111

-0.0889

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

-0.0882

0.3104

-0.0957

0.3313

-0.0472

n=

9

（5）取x(0)=（1,1，……，1,1）T，b=（1,1，……，1,1）T,迭代误差取默认值10 5，迭代次数上限取默认值50，使用Jacobi法进行迭代。

>>test2

>>b=zeros(20,1)

>>x0=zeros(20,1)

>>[x,n]=JacobiMethod(A,b,x0)

x=

0.2766

0.2327

0.2159

0.2223

0.2228

0.2221

0.2222

0.2222

0.2222

0.2222

0.2222

0.2222

0.2222

0.2222

0.2221

0.2228

0.2223

0.2159

0.2327

0.2766

n=

19

（6）取x(0)=（1,1，……，1,1）T，b=（1,1，……，1,1）T,迭代误差取默认值10 5，迭代次数上限取默认值50，使用Gauss-Seidel法进行迭代。

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

>>test2

>>b=zeros(20,1)

>>x0=zeros(20,1)

>>[x,n]=Gauss_SeidelMethod(A,b,x0)

x=

0.2766

0.2327

0.2159

0.2223

0.2228

0.2221

0.2222

0.2222

0.2222

0.2222

0.2222

0.2222

0.2222

0.2222

0.2221

0.2228

0.2223

0.2159

0.2327

0.2766

n=

10

(0)（7）取x=，b=Ae,将A的主对角线元素分别增大1,2,3……50倍，迭代误差取10 6，迭代次数上限取30，使用Jacobi法进行迭代，比较每次迭代次数。

>>n=test2_2()

n=

Columns1through30

21

6

55116595558585775656565656565

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

Columns31through50

5

5555555555555555555

2.2计算结果分析

由（1）~（6）中计算结果可见，一般情况下，使用Gauss-Seidel迭代法计算，其收敛速度要快于Jacobi迭代法。这是由于Gauss-Seidel迭代法在每次迭代过程中会使用最新的迭代值，因而收敛速度会比Jacobi迭代法更快。

由（7）中计算结果可见，随着不断增大A的主对角线元素，使用Jacobi迭代法计算的收敛迭代次数不断减小，直至减小到某一值才稳定下来。这是由于线性方程组的迭代收敛次数由B的范数决定，对于Jacobi法，B=D-1(L+U)，当A的主对角线元素不断增大时，B的范数不断减小，直至某一极限值，此时使用Jacobi迭代法计算的收敛迭代次数稳定到某一常数。

2.3M文件

2.3.1test2.m

%test2创建练习2中的系数矩阵A

clear;

e=ones(20,1);

A1=3*ones(1,20);

A2=-0.5*ones(1,19);

A3=-0.25*ones(1,18);

A=diag(A1)+diag(A2,1)+diag(A2,-1)+diag(A3,2)+diag(A3,-2);

.m2.3.2JacobiMethodJacobiMethod.m

function[X,n]=JacobiMethod(A,b,x0,eps,M)

%使用JacobiMethod函数求Ax=b的解

%A为系数矩阵

%b为常数向量

%x0为初值向量

%eps为给定误差界，可选

%M为限制迭代步数，可选

%X为线性方程组的解

%n为达到所需精度的迭代次数

%获取矩阵的阶数m

m=size(A);

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

m=m(1);

%对函数的输入参数进行判断处理

ifnargin==3

eps=1.0e-5;

M=50;

elseifnargin<3

error

return;

end

%判断A是否可逆，且对角元是否非0

flag=(m~=nnz(diag(A)));

if(det(A)==0)

disp('矩阵A不可逆，请重新输入');

error

return;

elseif(flag)

disp('矩阵A对角元有0元素，请重新输入');

error

return;

end

%求A的对角矩阵、下三角矩阵和上三角矩阵

D=diag((diag(A))');

L=tril(A)-D;

U=triu(A)-D;

%求等价方程组

Bj=D\(L+U);

f=D\b;

%创建元胞数组，用于存放迭代过程中的解向量

x=cell(1,M+1);

x{1}=x0;

x{2}=Bj*x{1}+f;

disp('x');

disp(1);

disp((x{1})');

disp('x');

disp(2);

disp((x{2})');

%迭代计算

fork=1:M

if(norm(x{k+1}-x{k},inf)>=eps)

x{k+2}=Bj*x{k+1}+f;

disp('x');

disp(k+2);

disp((x{k+2})');

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

elsebreak;

end

if(k+1==M)

disp('迭代次数超过给定值，可能不收敛！');

return;

end

end

%将迭代值赋给X,迭代次数赋给n

X=x{k+1};

n=k+1;

.m2.3.3Gauss_SeidelMethodGauss_SeidelMethod.m

function[X,n]=Gauss_SeidelMethod(A,b,x0,eps,M)

%使用Gauss_SeidelMethod函数求Ax=b的解

%A为系数矩阵

%b为常数向量

%x0为初值向量

%eps为给定误差界，可选

%M为限制迭代步数，可选

%X为线性方程组的解

%n为达到所需精度的迭代次数

%获取矩阵的阶数m

m=size(A);

m=m(1);

%对函数的输入参数进行判断处理

ifnargin==3

eps=1.0e-5;

M=50;

elseifnargin<3

error

return;

end

%判断A是否可逆，且对角元是否非0

flag=(m~=nnz(diag(A)));

if(det(A)==0)

disp('矩阵A不可逆，请重新输入');

error

return;

elseif(flag)

disp('矩阵A对角元有0元素，请重新输入');

error

return;

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

end

%求A的对角矩阵、下三角矩阵和上三角矩阵

D=diag((diag(A))');

L=tril(A)-D;

U=triu(A)-D;

%求等价方程组

Bg=(D-L)\U;

f=(D-L)\b;

%创建元胞数组，用于存放迭代过程中的解向量

x=cell(1,M+1);

x{1}=x0;

x{2}=Bg*x{1}+f;

disp('x');

disp(1);

disp((x{1})');

disp('x');

disp(2);

disp((x{2})');

%迭代计算

fork=1:M

if(norm(x{k+1}-x{k},inf)>=eps)

x{k+2}=Bg*x{k+1}+f;

disp('x');

disp(k+2);

disp((x{k+2})');

elsebreak;

end

if(k+1==M)

disp('迭代次数超过给定值，可能不收敛！');

return;

end

end

%将迭代值赋给X,迭代次数赋给n

X=x{k+1};

n=k+1;

2.3.4test2_2.mtest2_2.m

functionn=test2_2()

%用于求解习题2中第二问问题

%n是用来接收每次增长倍数相对应的迭代次数的向量组

clear;

e=ones(20,1);

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

A1=3*ones(1,20);

A2=-0.5*ones(1,19);

A3=-0.25*ones(1,18);

A=diag(A1)+diag(A2,1)+diag(A2,-1)+diag(A3,2)+diag(A3,-2);

b=A*e;

x0=zeros(20,1);

n=[];

fori=1:50

A1=i*3*ones(1,20);

A=diag(A1)+diag(A2,1)+diag(A2,-1)+diag(A3,2)+diag(A3,-2);

n(i)=Jacobi(A,b,x0,1e-6,30);

end

functionn=Jacobi(A,b,x0,eps,M)

%获取矩阵的阶数m

m=size(A);

m=m(1);

%对函数的输入参数进行判断处理

ifnargin==3

eps=1.0e-5;

M=50;

elseifnargin<3

error

return;

end

%判断A是否可逆，且对角元是否非0

flag=(m~=nnz(diag(A)));

if(det(A)==0)

disp('矩阵A不可逆，请重新输入');

error

return;

elseif(flag)

disp('矩阵A对角元有0元素，请重新输入');

error

return;

end

%求A的对角矩阵、下三角矩阵和上三角矩阵

D=diag((diag(A))');

L=tril(A)-D;

U=triu(A)-D;

%求等价方程组

Bj=D\(L+U);

f=D\b;

%创建元胞数组，用于存放迭代过程中的解向量

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

x=cell(1,M+1);

x{1}=x0;

x{2}=Bj*x{1}+f;

%迭代计算

fork=1:M

if(norm(x{k+1}-x{k},inf)>=eps)

x{k+2}=Bj*x{k+1}+f;

elsebreak;

end

if(k+1==M)

disp('迭代次数超过给定值，可能不收敛！');

return;

end

end

%将迭代值赋给X,迭代次数赋给n

X=x{k+1};

n=k+1;

习题3

用迭代法求方程x+3x 1=0的全部根，要求误差限为0.5×10。32 8

3.1程序

（1）通过作图找到有根区间为（-3，-2.5），（-1，-0.5),(0.5,1)

>>test3

>>

plot(x,y)

图1非线性方程组单根区间

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

（2）使用二分法求方程根

>>[root,n]=halfIntervalMethod(f,-3,-2.5,0.5e-8)

root=

-2.8794

n=

27

>>[root,n]=halfIntervalMethod(f,-1,-0.5,0.5e-8)

root=

-0.6527

n=

27

>>[root,n]=halfIntervalMethod(f,-0.5,1,0.5e-8)

root=

0.5321

n=

29

（3）先使用二分法求方程近似根，再将其作为初值，使用牛顿法求解（此处只求区间（-3，-2,5）上的根）

>>[root,n]=halfIntervalMethod(f,-3,-2.5,1e-2)

root=

-2.8789

n=

6

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

>>[root,n]=NewtonMethod(f,-2.8789,0.5e-8,10)

k=

1

-2.8794

-2.8794

root=

-2.8794

n=

3

3.2计算结果分析

由计算结果可见，只使用二分法求区间（-3，-2,5）上的根，达到误差限要求时所进行的迭代步数为27，而使用二分法与牛顿法的组合则一共仅用9步，这是因为牛顿法至少是2阶收敛，显然效率更高。

3.3M文件

3.3.1test3.m

%test3用于创建习题3中运算表达式

clear;

symsx

f=x^3+3*x^2-1;

%通过作图法找到有根区间为（-3，-2.5），（-1，-0.5),(0.5,1)

x=-3.5:0.1:1;

y=subs(f,x);

.m3.3.2halfIntervalMethodhalfIntervalMethod.m

function[root,n]=halfIntervalMethod(f,a,b,eps)

%二分法求非线性方程f（x）=0在某个单根区间上的一个零点

%f为非线性方程的函数表达式

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

%a，b分别为所选求根区间的左右端点

%eps为给定误差界，可选

%root为所求区间上的单根

%n为迭代次数

%对输入参数进行判断处理

if(nargin==3)

eps=1.0e-3;

end

%赋初值

fa=subs(f,a);

fb=subs(f,b);

root=(a+b)/2;

fr=subs(f,root);

n=0;

%用二分法求根

while(abs(a-b)>eps)

if(fa*fr<0)

b=root;

root=(a+b)/2;

elseif(fb*fr<0)

a=root;

root=(a+b)/2;

end

fa=subs(f,a);

fb=subs(f,b);

fr=subs(f,root);

n=n+1;

end

.m3.3.3NewtonMethodNewtonMethod.m

function[root,n]=NewtonMethod(f,x0,eps,M)

%NewtonMethod用牛顿法，即切线法求非线性方程f（x）=0在某个单根区间上的一个零点

%f为非线性方程的函数表达式

%x0

%eps为给定误差界，可选

%M为迭代次数限定值

%root为所求区间上的单根

%n为迭代次数

%对输入参数进行判断处理

if(nargin==3)

本报告内容为大连理工大学研究生院2012年矩阵与数值分析课程上机作业代码及报告分析,内容丰富详实,且全部由作者独立完成,没有任何剽窃,程序由Matlab编写完成,希望对大家有所帮助!

eps=1e-5;

M=50;

end

%计算参数初始化

dfdx=diff(f);

x=[];

x(1)=x0;

x(2)=x(1)-subs(f,x(1))/subs(dfdx,x(1));

k=1

%迭代求根

while(abs(x(k+1)-x(k))>eps)

x(k+2)=x(k+1)-subs(f,x(k+1))/subs(dfdx,x(k+1));

disp(x(k+1));

if(k+2>M)

disp('迭代次数超过限定值，迭代过程可能不收敛');

return

end

k=k+1;

end

root=x(k+1);

n=k;

习题4

给定数据如下表：

xi

yi00.010.7921.5332.1942.7153.0363.2772.8983.0693.19103.29编制程序求三次样条插值函数在插值中点的样条函数插值，并作点集{xi,yi}和样条插值函数的图形，满足的第一、二边界条件为

s′(0)=0.8,s′(10)=0.2.s′′(0)=s′′(10)=0.

4.1程序

（1）利用第一边界条件求数据点集的三次样条插值函数

>>test4

>>s=tripleSpindleMethod(X,Y,0.8,0.2,1)

本文来源：https://www.bwwdw.com/article/92ae.html