数值分析上机作业

更新时间：2023-07-17 11:09:01 阅读量：实用文档文档下载

说明：文章内容仅供预览，部分内容可能不全。下载后的文档，内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的，是否完整无缺。

数值是什么意思推荐度：
相关推荐

姓名：学号：学院：专业：

《数值分析》上机作业

4题

设计程序如下： clear;

%定义x，y的区间

x=-10:0.2:10; %步长为100 y=-10:0.2:10;

%在XOY平面生成二维网格数据 [X,Y]=meshgrid(x,y); %对二元函数进行表达

a=-abs(X);b=X+Y;c=X.^2+Y.^2+1; Z=exp(a)+cos(b)+1./c; %绘制三维图形 mesh(X,Y

,Z);

图1. 区间等分100份图2.区间等分200份

图3.区间等分400份

2题

设计程序如下： Clear

n=8; V=220; R=27;

a=[0 -2 -2 -2 -2 -2 -2 -2]; b=[2 5 5 5 5 5 5 5]; c=[-2 -2 -2 -2 -2 -2 -2]; I=[V/R 0 0 0 0 0 0 0]; for i=2:n

a(i)=a(i)/b(i-1); b(i)=b(i)-c(i-1)*a(i); I(i)=I(i)-a(i)*I(i-1); end;

I(n)=I(n)/b(n); for i=n-1:-1:1

I(i)=(I(i)-c(i)*I(i+1))/b(i); end; I;

fprintf(“各电路的电流量I为”); 运行结果如下：

各电路的电流量I为

I=8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194

0.0477

2题

高斯赛德尔迭代

设计程序如下：

：A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;] x=[0;0;0;0;0]; b=[12;-27;14;-17;12] c=0.000001 L=-tril(A,-1)

U=-triu(A,1) D=(diag(diag(A)))

X=inv(D-L)*U*x+inv(D-L)*b; k=1;

while norm(X-x,inf)>= c x=X;

X=inv(D-L)*U*x+inv(D-L)*b; k=k+1; end X k

计算结果：X = 1.0000 -2.0000 3.0000 -2.0000 1.0000 k =37

Jacobi迭代法：设计程序如下：

Jacobi：b=[12;-27;14;-17;12] x = [0;0;0;0;0;] k = 0; r = 1; e=0.000001

A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;] D = diag(diag(A)); B = inv(D)*(D-A); f = inv(D)*b;

p = max(abs(eig(B))); if p >= 1

'迭代法不收敛' return end

while r >e x0 = x;

x = B*x0 + f; k = k + 1;

r = norm (x-x0,inf); end x k

计算结果：x = 1.0000 -2.0000 3.0000 -2.0000 1.0000 k =65

SOR：A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15] x=[0;0;0;0;0]; b=[12;-27;14;-17;12] e=0.000001

w=1.44; L=-tril(A,-1)

U=-triu(A,1) D=(diag(diag(A)))

X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*b n=1;

while norm(X-x,inf)>=e x=X;

X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*b; n=n+1; end X n

计算结果：X = 1.0000 -2.0000 3.0000 -2.0000 1.0000 n =22

结果分析：由迭代次数的比较可知该情况下Jacobi迭代法比Gauss-Seidel迭代法收敛的慢

2题

设计程序如下

%幂法和反幂法

A=[12 6 -6;6 16 2;-6 2 16]; x0=[1;1;1];

%% 幂法求按模最大特征值 x=x0;

N = 3;

lamada = max(x0); for k =1:N

y = x/lamada; x = A*y;

lamada = max(x); end lamada

%% 反幂法求精确值 lamda0 = floor(lamada); I = [1 0 0;0 1 0;0 0 1]; [L,U,P] = lu(A-lamda0*I);

x=x0; k=0; error = 1;

maxe = 1e-10; alpha = max(x0); beta = alpha;

while error>maxe alpha = beta; y = x/alpha; z = L\y; x = U\z;

beta = max(abs(x));

error = abs(1/beta-1/alpha); k=k+1; if(k==100) break; end end

lamda = lamda0+1/beta k

1题

设计程序如下：

subplot(2,1,1) syms x

for i=1:10

for k=11:-1:i+1 f(k)=(f(k)-f(k-1))/(i); end end

a=-5:0.1:5;

b=1./(1+4.*(a).^2); plot(a,b,'+');

title('\fontname{隶书} 原函数和插值多项式') hold on

p=0; for i=1:11 q=1;

for j=1:1:i-1 q=q*(x-j+6); end; p=p+f(i)*q; end; p

axis([-5,5,-1,5]);

x=-5:0.1:5; e=subs(p) plot(x,e);

hold on

text(2,3,'+++插值函数') text(2,2.2,'-----原函数')

subplot(2,1,2)

for i=-5:5

f(i+6)=1/(1+4*(i)^2); end

e-b

plot(x,e-b); hold on;

title('\fontname{隶书} 误差图') p =

(36*x)/6565

(3550298616520539*(x

4)*(x

5))/1152921504606846976

(2689247898264063*(x + 3)*(x + 4)*(x + 5))/1152921504606846976 + (1806978031308661*(x + 2)*(x + 3)*(x + 4)*(x + 5))/576460752303423488 + (462354082176629*(x + 1)*(x + 2)*(x + 3)*(x + 4)*(x + 5))/144115188075855872 - (5850230976024283*x*(x + 1)*(x + 2)*(x + 3)*(x + 4)*(x + 5))/1152921504606846976 + (6518522480310501*x*(x - 1)*(x + 1)*(x + 2)*(x + 3)*(x + 4)*(x + 5))/2305843009213693952 - (2258610859405831*x*(x - 1)*(x + 1)*(x - 2)*(x + 2)*(x + 3)*(x + 4)*(x + 5))/2305843009213693952 + (1143600435142193*x*(x - 1)*(x + 1)*(x - 2)*(x + 2)*(x - 3)*(x + 3)*(x + 4)*(x + 5))/4611686018427387904 - (7319042784910035*x*(x - 1)*(x + 1)*(x - 2)*(x + 2)*(x - 3)*(x + 3)*(x - 4)*(x + 4)*(x + 5))/147573952589676412928 + 49/1313

2题

设计的程序如下： Clear

X=[2 3 5 6 7 9 10 11 12 14 16 17 19 20];

Y=[106.42 108.26 109.58 109.50 109.86 110.00 109.93 110.59 110.60 110.72 110.90 110.76 111.10 111.30]; D=1./X; S=1./Y; S=S';

A=ones(14,2); for i=1:14

A(i,2)=D(i); end; c=A\S; a=c(1); b=c(2); syms x;

y=x/(a*x+b);

disp('y=x/(a*x+b)'); a

b x1=[1:0.1:21]; y1=x1./(a*x1+b);

plot(X,Y,'*',x1,y1,'k'); title('最小二乘拟合函数');

legend(' 原始点','y=x/(a*x+b)');

运行结果如下： y=x/(a*x+b) a =0.0090

b = 8.4169e-004

1题

设计的程序如下：

a=0; X=[]; for b=-5:0.05:5 w=1.0e-3;

F1=cos(a)+cos(b^2/2); F2=cos((a+b)^2/8);

S0=((b-a)/6)*(F1+4*F2); m=2; h=(b-a)/4; F3=0;

for k=1:2^(m-1)

F3=F3+cos((a+(2*k-1)*h)^2/2); end

S=(h/3)*(F1+2*F2+4*F3); while abs(S-S0)>15*w m=m+1;

h=h/2; F2=F2+F3; S0=S; F3=0;

for k=1:2^(m-1)

F3=F3+cos((a+(2*k-1)*h)^2/2); end

S=(h/3)*(F1+2*F2+4*F3); end

X=[X S]; end

Y=[]; for b=-5:0.05:5

w=1.0e-7;

F1=sin(a)+sin(b^2/2); F2=sin((a+b)^2/2); S0=((b-a)/6)*(F1+4*F2); m=2; h=(b-a)/4; F3=0;

for k=1:2^(m-1)

F3=F3+sin(a+(2*k-1)*h)^2/2; end

S=(h/3)*(F1+2*F2+4*F3); while abs(S-S0)>15*w m=m+1; h=h/2; F2=F2+F3; S0=S; F3=0;

for k=1:2^(m-1)

F3=F3+sin((a+(2*k-1)*h)^2/2); end

S=(h/3)*(F1+2*F2+4*F3); end

Y=[Y S]; end

plot(X,Y

,'r.-')

第八章

设计的程序如下：

a=0.5;

b=128.52*log((513+0.6651*a)/(513-0.6651*a)); while(abs(a-b)>1.0e-8) a=b;

b=128.52*log((513+0.6651*a)/(513-0.6651*a)); end

第九章

1题

Adams预测-校正设计的程序如下：

a=0; b=3.14; N=30; y(1)=1; h=(b-a)/3;

x(1)=a; n=2;

while n<5

f(n-1)=-y(n-1)+2*cos(x(n-1));

K1=h*f(n-1);

K2=h*(-(y(n-1)+K1/2)+2*cos(x(n-1)+h/2)); K3=h*(-(y(n-1)+K2/2)+2*cos(x(n-1)+h/2)); K4=h*(-(y(n-1)+K3)+2*cos(x(n-1)+h)); y(n)=y(n-1)+(1/6)*(K1+2*K2+2*K3+K4); x(n)=a+(n-1)*h;

f(n)=-y(n)+2*cos(x(n)); n=n+1; end p(4)=0; c(4)=0; while n<30

x(n)=x(n-1)+h;

f(n-1)=-y(n-1)+2*cos(x(n-1));

p(n)=y(n-1)+h/24*(55*f(n-1)-59*f(n-2)+37*f(n-3)-9*f(n-4)); m(n)=p(n)+251/270*(c(n-1)-p(n-1));

c(n)=y(n-1)+h/24*(9*(-m(n)+2*cos(x(n)))+19*f(n-1)-5*f(n-2)+f(n-3)); y(n)=c(n)-19/270*(c(n)-p(n)); n=n+1;

end

for i=1:1:29

err(i)=sin(x(i))+cos(x(i))-y(i); end x,y,err x =

Columns 1 through 11

0 1.0467 2.0933 3.1400 4.1867 5.2333 6.2800 8.3733 9.4200 10.4667

Columns 12 through 22

11.5133 12.5600 13.6067 14.6533 15.7000 16.7467 17.7933 19.8867 20.9333 21.9800