计算方法上机实验报告-MATLAB

  • 格式:docx
  • 大小:721.02 KB
  • 文档页数:18

下载文档原格式

  / 18
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

《计算方法》实验报告

指导教师:

学院:

班级:

团队成员:

一、题目

例2.7应用Newton 迭代法求方程210x x --=在1x =附近的数值解

k x ,并使其满足8110k k x x ---<

原理:

在方程()0f x =解的隔离区间[],a b 上选取合适的迭代初值0x ,过曲线()y f x =的点()()

00x f x ,引切线

()()()1000:'l y f x f x x x =+-

其与x 轴相交于点:()()

0100 'f x x x f x =-,进一步,过曲线()y f x =的

点()()11x f x , 引切线

()()()2111: 'l y f x f x x x =+-

其与x 轴相交于点:()

()

1211 'f x x x f x =-

如此循环往复,可得一列逼近方程()0f x =精确解*

x 的点

01k x x x ,,,,,其一般表达式为:

()()

111 'k k k k f x x x f x ---=-

该公式所表述的求解方法称为Newton 迭代法或切线法。

程序:

function y=f(x)%定义原函数

y=x^3-x-1;

end

function y1=f1(x0)%求导函数在x0点的值

syms x;

t=diff(f(x),x);

y1=subs(t,x,x0);

end

function newton_iteration(x0,tol)%输入初始迭代点x0及精度tol x1=x0-f(x0)/f1(x0);k=1;%调用f函数和f1函数

while abs(x1-x0)>=tol

x0=x1;x1=x0-f(x0)/f1(x0);k=k+1;

end

fprintf('满足精度要求的数值为x(%d)=%1.16g\n',k,x1);

fprintf('迭代次数为k=%d\n',k);

end

结果:

二、题目

例3.7试利用Jacobi 迭代公式求解方程组

1234451111101112115181

11

1034x x x x ----⎛⎫⎛⎫

⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪= ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪---⎝⎭⎝⎭

⎝⎭ 要求数值解()k X 满足()42

10k --≤X X ,其中T (1,2,3,4)=X 为方程

组的精确解。

原理:

将线性方程组的系数矩阵1112121

2221

2

n n n n nn a a a a a a a a a ⎛⎫ ⎪

⎪= ⎪ ⎪⎝⎭

A 分解为=++A L D U ,其中1122=(,,,)nn diag a a a D ,

31

32,1

211

2

000000000=0n n n n a a a a a a -⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭

L ,12

1312321,00=00

00000

0n n n n a a a a a a -⎛⎫ ⎪

⎪ ⎪ ⎪ ⎪⎝⎭

U 当对角阵D 可逆时,方程组=AX b 可等价地写成

()11 --=-++X D L U X D b ,

据此可得Jacobi 迭代公式为

(

)

()()111 0,1,

k k

k +--=-++=X D L U X D b, ,

其中()

()()()

(

)

T

12=,k k k n n

k x x x ∈X

R ,,,该迭代公式也可以写成如下的分

量形式

1(()11)

,,1,2,,1n

k k i ij j j j i ii i

x

b a i n a x ++=≠⎛⎫= ⎪⎭

=⎝-∑

程序:

function jacobi()

%输入矩阵A 、b 、精度tol%

A=[5 -1 -1 -1; -1 10 -1 -1; -1 -1 5 -1; -1 -1 -1 10]; b=[-4 12 8 34]'; tol=10^-4;

% A=input('系数矩阵A=');%输入矩阵A % b=input('矩阵b= ');%输入矩阵b

% tol=input('精度要求tol=');%输入精度tol X=inv(A)*b; [n,~]=size(A); D=diag(diag(A)); L=tril(A)-D; U=triu(A)-D; X0=zeros(n,1);

X1=inv(D)*b-inv(D)*(L+U)*X0; X0=X1; X2=X-X0; k=1;

while abs(norm(X2,2))>=tol

X1=inv(D)*b-inv(D)*(L+U)*X0; X0=X1; X3=X-X0; %收敛性判断%

if norm(X3,2)>norm(X2,2) break; else

X2=X3; k=k+1; end end

% 结果输出%

if X2~=X3