华科计算方法上机作业代码

《计算方法》上机实验试题

1. (25分)计算积分

dx x x I n

n ⎰

+=1

5

， n=0,1,2,…,20 若用下列两种算法 (A) n

I I n n 151+-=- (B) ⎪⎭

⎫ ⎝⎛-=

-n n I n I 1511 试依据积分I n 的性质及数值结果说明何种算法更合理。

i0=1.8232000e-001; disp(i0); for n=1:30

i1=double(-5.0*i0+1.0/n); disp(n); disp(i1); i0=i1; end

0.1823

1

0.0884

2

0.0580

3

0.0433

4

0.0333

5

0.0333

6

2.4261e-13 7

0.1429

8

-0.5893

9

3.0575

10

-15.1877

11

76.0294

12

-380.0637

13

1.9004e+03 14

-9.5019e+03 15

4.7510e+04 16

-2.3755e+05 17

1.1877e+06 18

-5.9387e+06 19

2.9693e+07 20

-1.4847e+08 21

7.4234e+08 22

-3.7117e+09 23

1.8558e+10 24

-9.2792e+10 25

4.6396e+11 26

-2.3198e+12 27

1.1599e+13 28

-5.7995e+13 29

2.8998e+14 30

-1.4499e+15

i30=5.9140e-003;

disp(i30);

for n=30:-1:1

i1=double(-0.2*i30+1.0/(5*n)); disp(n);

disp(i1);

i30=i1;

end

0.0059

30

0.0055

29

0.0058

28

0.0060

27

0.0062

26

0.0065

25

0.0067

24

0.0070

23

0.0073

22

21 0.0080 20 0.0084 19 0.0088 18 0.0093 17 0.0099 16 0.0105 15 0.0112 14 0.0120 13 0.0130 12 0.0141 11

10

0.0169

9

0.0188

8

0.0212

7

0.0243

6

0.0285

5

0.0343

4

0.0431

3

0.0580

2

0.0884

1

0.1823

2.(25分)求解方程f(x)=0有如下牛顿迭代公式

)

()

(111---'-

=n n n n x f x f x x ， n ≥1，x 0给定

(1)

编制上述算法的通用程序，并以ε<--1n n x x （ε为预定精度）作为终止迭代的准则。

(2)

利用上述算法求解方程

0cos :)(=-=x x x f

这里给定x 0=π/4，且预定精度ε=10-10。

fun=inline('cos(x)-x')

fun =

Inline function: fun(x) = cos(x)-x

>> dfun=inline('-sin(x)-1')

dfun =

Inline function: dfun(x) = -sin(x)-1

>> x=agui_newton(fun,dfun,pi/4,1e-10) 0.7395

0.7391

0.7391

0.7391 x =

0.7391

3. (25分) 已知)3()(x x e x e x f -=，

(1)

利用插值节点x 0=1.00，x 1=1.02，x 2=1.04，x 3=1.06，构造三次Lagrange 插值公式，由此计算f(1.03)的近似值，并给出其实际误差； function f = Language(x,y,x0) syms t;

if(length(x) == length(y)) n = length(x); else

disp('x 和y 的维数不相等！'); return;

end %检错

f = 0.0; for(i = 1:n) l = y(i); for(j = 1:i-1)

l = l*(t-x(j))/(x(i)-x(j)); end;

for(j = i+1:n)

l = l*(t-x(j))/(x(i)-x(j)); %计算拉格朗日基函数 end;

f = f + l; %计算拉格朗日插值函数 simplify(f); %化简

if(i==n)

if(nargin == 3)

f = subs(f,'t',x0); %计算插值点的函数值 else

f = collect(f); %将插值多项式展开

f = vpa(f,6); %将插值多项式的系数化成6位精度的小数 end end end

function f=fun(x) f=exp(x)*(3*x-exp(x)); end

x0=[1.00 1.02 1.04 1.06] x0 =