数值分析实验报告一

实验报告

课程名称：_________ 数值分析________________ 实验项目：解方程/解方程组

专业班级：_________ 机械工程_________________ 姓名：李奇学号:201820065

实验室号：___________ 实验组号：___________ 实验时间：___________ 批阅时间：___________ 指导教师: 曲绍波成绩: _______________

沈阳工业大学实验报告

(适用计算机程序设计类)

专业班级：机械工程学号：201820065 姓名：李奇实验名称：解方程

沈阳工业大学实验报告

附件A

（适用计算机程序设计类）

专业班级：机械工程学号：201820065 姓名：李奇实验步骤或程序：

1、编写牛顿法M文件

山扁辑諾 * C:\U s ers\l i qi\De s ktop\n ewton, m (?) newton￥m +

1 - syms f x:

2 -f= input（请输入）；

3 -df^diff(f);

4 _K0-input ＜请输入迭代初值;

5 -e= input （' iWSi入根的误差限亡才）；

6 -input C请輸入迭代次数限N="）;

1 N k=l ;

S - while (k

9 -J5=^ 0 ；

10 ->:L=xO-eval (f )/eval (df):

11 -if abs(Kl-5￡0)

12 -fprintf t r沪k 3f\ri迭f?次数为r xl.k)

13 -break

14 -else

15 -K0=K1;

16 -k=k+l;

17 - 巳nd

15 -end

19 -if k>=N

20 - fprintf (失B&\n')

21 一end

2、调用newton.m文件，输入求解

带令行圖口

? neuron

fx请输入f

3、运行结果如下

?? neirt on

请输入f ) =cos (K l-z*exp (x)

i青输入迭代初值菠oo

语输入根的误姜限e=10^5

请输入迭优次数限

x-0.517757

迭代衣數为用

4、编写割线法M文件

(1子程序

N 秦蛊藹-C:\User5\liq i\De5ktDp\fun. m

fun.n

I 2

1C '

11 12 13function K= secant(xO,xl T tol)

if nar gin < 10 tol N 1. Oe-6;

x = xO - fun(xO) * (xO - xl) / (fun(xO) - fun(xl)):

n = 1:

我订已(abs(xO-xl) > tol)蜒(n <- 1000)

xl = xO;

xO = x:

x = x0 - fun@Q) * (xO - xl) / (fun(xO) - fun(xl)); n = n + 1; end

5、调用secant.m文件，输入方程

命令行窗口

>> secant(0,0. 2*le-6)

&运行结果如下

命令行窗口

>> secant (Os 0. 2, le~6)

n =

8

arts =

0.5178

沈阳工业大学实验报告

（适用计算机程序设计类）

专业班级：机械工程学号：201820065 姓名:李奇

实验名称：解方程组

1. 实验目的:

掌握高斯消去法，列主元高斯消去法和 LU 分解法求解方程组 能够使用MATLA 编程

2. 实验内容:

（1）学习掌握高斯消元，列主元，LU 分解原理 ⑵ 学习MATLA 语言的语法格式。 （3）编写MATLA 程序求解下列方程组 4x 5.3y -5.6z -3m -3.4n =100.16 5x-2.1y 3.2z 4m -8n =-75.72

5x-3y -8z 2.3m 3n =57.1 4.2x-3y -2n =3.72

3. 实验方案（程序设计说明）

高斯消去法的基本思想是通过消元把线性方程组化为等价的上三角方程组，再 进行求解。高斯消去一般由“消元过程”和“回代过程“两部分组成。消元过程 就是按确定的计算过程对方程组的增广矩阵进行初等变换，将原方程化为与之等 价的上三角方程组；回代过程就是对得到的上三角方程组求解的过程。

列主元高斯消去法的算法步骤：（1）对k=1,2,…,n -1，选主元素，选出｛|a （k ）

kk

|，|a （k

）k+1，k |,…，|a （k ）nk | ｝中绝对值最大的兀素|a （k ）mk |，将第k 行和第m 行交换。（2）

消 元过程。（3）回代过程。⑷ 输出方程组的解。

LU 分解法的算法步骤：（1）首先进行LU 分解，计算U 的第一行和L 的第一列。 ⑵对k=2, 3,…，n,计算U 的第k 行和L 的第k 列。⑷计算y,求解Ly=b 。 ⑷求解x,求解Ux=y 。

4. 实验步骤或程序（经调试后正确的源程序）

见附件B

5. 程序运行结果

见附件B

6. 出现的问题及解决方法

II 2x-4y-7.2z-5m-2.4n =98.2

实验步骤或程序:

1、编写高斯消元法M 文件

-C:\U sers\liqi\D?Ektop\Gaus5xi3o qu.m Gaussxiaoqu.rn +

1

. tunct ion 2RA h n. Xj -jausssiaoqutAj bj

2 - S=：A b];

3 - n=length

4 - RA=rank (A):

RE=rank'.Ej :

if RB-RA 汕

di E pC^示：因均肥、眼 蘭以it 右程坦无詮’i

r et urr

if RA==RE

i￡ R^-n

dispC4S^ ：因為RpRBgfif 以哉方程组有唯一的解’) X^zeros (n. 1): fcr p= 1: n-1

.-

f ar k=p-i-l: r

JD =B (k, p>/'E (p ： p) : B (k, p: n+I ) -B (k> p. r.+ l)-m*Etp n p:

n+1): end ^nd

W(L n T n+l):A-E(l;n T l:n) ;X

for

. _

1. 1

X t a) = (b (o)-siin (A(u. q+1 n) *X (Q +1 : n) 1 )/A(Q , Q )

€nd else

cisp(^提示：因为fUB 恥5,所以此方程通自无穷毋蔚’)

end

end

2、调用Gaussxiaoqu.m 文件，输入方程组

盛Ci 雲匚

?

>>

I aaqu ?臼丨

R -3.^ ^-2.

隊7 -&2.1 3:4. 2.-3.D. 0.^23:

? b=\O0,H -"C. 72^5,2 -67.1;3,7?'.: fit ? &刊垢3>： z 四m

附件B

专业班级： 机械工程

沈阳工业大学实验报告

（适用计算机程序设计类）

学号： 201820065

姓名： 李奇

Z 廉数担弾，S 酣扫端项

牺拘摺广担痔

t 右端项b 的维数

妇系數矩降呦秩 比培广拒蹲吕审秩

10 -

11 -

12 - 13 - 14 -

!□-

1￡ - 17 - 18 - 19 - 20 -

St -

22 - 23 - 24 - 25 - 26 -

6 -

7 - S -

3、运行结果如下

>> taussxiaoqu bj A

提示：因为RA=RE=n）所以.此方程组有唯一的解

-11. 0000

3.0000

-7. 5000

-lk 0000

3. 0000

(|

3.2000 -7,5000 -11.

0000

3.0000

1. 6000

3, 2000

-7.6000

-11. 0000

3, OOOG

ans =

5、调用Gauss.m 文件，输入方程组

齒甘右誉口

?

?如」期心为)

>? A-：l, 6, ~3> -3. 4; 5.-2. X 3■签 4-$;岛一￡ 一訴監-乩 TY 壬一為 f Z. 3, 3 h J K 2.-J. C* C. -2：

…14 -T5. 72 ?.f;97.1 ;3. 72]

Jx Gauss ti 1

4、编写高斯列主元法M 文件

函

- C:\Userj\liqi\0e5 ktop\Ga uss

/n

G RU 販m

+

1

Fl function LKA.BE, r.?X?=>ajss (A, t)

2 - B- Ci b]:

3 - n 2

-length(b i ; 4 -

RA-r^nk (A): %A 为爲新靈偉，b 为右旅质—

弔E 拘増广矩降 务右端项b 用维藪 舌辛柚垢畛的徘

3 - if cha>0

9 e dispC 提示：因为肥一我.所以此方谨俎无解「) 10 r r^t nrn 11 a end

12 -

if R.fl

==RB

13 - 14 -

坊-

1：-

13 -

19 -

i ￡ RA.??JI

disp('^^；匡畑=M=n.所以此方程组苗哇一胡硒' ) X=3eios (r T 1! : ?=ZSTOE (l n rrHl 1

; for p_i ； n~ 1

[y, j]-toM(abs Etp:n t p)i) 勺戎出歹」屮绽刈直最玄

的對 C-E (p,:):

20

21

22

刀

24 25

商

2'

23

23 3(1

3L

32

33

B | j*p'l s :)-C;

fcr k-p^-1: u 乜肖兀

nr-E(k f p)/E (p E p) , B Oc f p ntl )-B(k f p: n+1) ■-s*B (p t p:irtl), ftr.d

end

b=B (1:11, n+1 J ; JkB (1: n, 1 :n ::(nJ /A (a. n) :

5回仕

for c^n-1 ： -1: 1

X(4)- (b (q) -sum 4A q t q-l-1: n) *X(Q ^1 : B ) 1

g);;

end

end #ls?

dLsp （提示：S.X]RA=FE

ei:d

&运行结果如下

-11. -1(1

fCJ

2. OOCO

-弓0〔〔1

-11. 30C0

3. jflCO

3. 2Q

〔〔I

-7. 5QCO

11. oooo

3. DOCO

1.6Q00

3.2000

5000

-JEOOOO

3,0000 Ans

7、编写LU分解法M文件

N 索垠辭-C:\ll s ers\l iq i\De sktc p \ L U. im

LU, r n X 1 +

1E■] function LU(A? bi—I 2.n. Ji.-S12€ (Al :%靭始比矩陆丄?HL和U

3n= length(h);

4L= aye (n)1:酚前僮陆

5U=2?ros(n):为n翫番方時

c if n^-m

$di3p(7

"不等于”无法逛行LU分琳"i

8

return;

9ct;d

10-fcx . n-l弔判断!广1时顺序主子戒是否

対零

11if (det (A(l：i, 1: i) )==0)

12disp(^提冷矩降阳f在为零的噸予圭于式,无法进行LU讣蔚)

13return;

14else

15

end

17end

18

J f cr j-1: n hU的第一行■郵第一行19U(L j)-A(U)；

20er,d

21 1 fcr i=2; n紅的第一列屯的第一列22Ki, 1)/UH T 1);

23end

24fk-2:n

25

U

26

L 0+1: iuk)-(A(k+1 ：n,k)-L(k+l：na：l^l *U(1 k-l7kJ)/U (k T k);谕>」尊岀L 的第27end

2S L

29U需铺出1：矩阵

30g 1);片莎皓褂方程齟Ur-y 31;

32—￡ i>i k-2. n

33_ 1 k.1-1. 'k 1 - bum L 1 k- 1 1 +_- 1 k-l'1■

34end

3o y

36x

37-El fir k-n-l:-'l： I

3S x(k)= iy

39end

4G X

8、调用LU.m 文件，输入方程组

LJ (A, bi

9、运行结果如下

1.

0000 0 0

0 0 E2500 LOOOO 0

0 0.5000 0.7522 EOOOO 0

0 k2500

1.1032 1. 0064 1. OOGO 0 1.0500

0.9317

0. 3395

-0.1815

1. 0000

u =

4.0000

5.3000 -乱 6000

-3. 0000 -3. 4000 0 ^-8.7250

10. 2000 7,7500 -3… 7500 0 0 -12. 1742

-9.4069 2.1582 0 0 0 6. 9677

9. 2148 0

6.1907

IDO.1600 -200.9200 20E2567 -48. &999 18.5722

5o2.0549 206.8016 24. 6226 -31. 5945 3. 0000

L

＞ 4, S. 3. ~5? 6, _3i-3. 4■: &b _

2,1, 3. & 4. ~￡ ； 2,

? BS K 1B;-75*72;S8.2;57.1:3.72：;

/i

、＞ LU(A,b)

一＞k 一 ；■ 2* -乩一2? 4 ； 5i -乳 一 & 2. J. M 4?畫一丄(L th -2】