数值分析上机题目4

数值分析上机题目4 -CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN

实验一

实验项目：共轭梯度法求解对称正定的线性方程组 实验内容：用共轭梯度法求解下面方程组

(1) 123421003131020141100155x x x x -⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪=

⎪ ⎪ ⎪-- ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ 迭代20次或满足()(1)

1110k k x x --∞

-<时停止计算。

编制程序：储存m 文件

function [x,k]=CGmethod(A,b)

n=length(A);x=2*ones(n,1);r=b-A*x;rho=r'*r; k=0;

while rho>10^(-11) & k<1000 k=k+1; if k==1 p=r; else

beta=rho/rho1; p=r+beta*p; end w=A*p;

alpha=rho/(p'*w); x=x+alpha*p; r=r-alpha*w; rho1=rho; rho=r'*r; end

运行程序：

clear,clc

A=[2 -1 0 0;-1 3 -1 0;0 -1 4 -1;0 0 -1 5]; b=[3 -2 1 5]'; [x,k]=CGmethod(A,b)

运行结果： x =

(2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵， A[i,i]=4，A[i,i-1]=A[i-1,i]=-1，i=2,3,..,n b[1]=3, b[n]=3, b[i]=2，i=2,3,…,n-1

迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。 编制程序：储存m 文件

function [x,k]=CGmethod_1(A,b) n=length(A);x(1:n,1)=0;r=b-A*x;r1=r; k=0;

while norm(r1,1)>10^(-7)&k<10^4 k=k+1; if k==1 p=r; else

beta=(r1'*r1)/(r'*r);p=r1+beta*p; end r=r1; w=A*p;

alpha=(r'*r)/(p'*w); x=x+alpha*p; r1=r-alpha*w; end

运行程序： clear,clc n=1000; A=hilb(n); b=sum(A')';

[x,k]=CGmethod(A,b)

实验二

1、 实验目的：用复化Simpson 方法、自适应复化梯形方法和Romberg 方法求数值积分。

实验内容：计算下列定积分

(1) dx x x x ⎰⎪⎪⎭

⎫ ⎝⎛+-2

02

610 (2)⎰10dx x x (3) ⎰20051dx x 实验要求：

（1）分别用复化Simpson 公式、自适应复化梯形公式计算要求绝对误差限为7102

1

-⨯=

ε，输出每种方法所需的节点数和积分近似值，对于自适应方法，显示实际计算节点上离散函数值的分布图；

（2）分析比较计算结果。 程序： syms x

f=x^6/10-x^2+x %定义函数f （x ） n=input('输入所求导数阶数:')

f2=diff(f,x,n) %求f(x)的n 阶导数 （1）复化梯形 clc clear

syms x%定义自变量x

f=inline('x^6/10-x^2+x','x') %定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可

f2=inline('3*x^4 - 2','x') %定义f(x)的二阶导数,输入程序1里求出的f2即可。

f3='-(3*x^4 - 2)'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值e=*10^(-7) %精度要求值

a=0 %积分下限

b=2 %积分上限

x1=fminbnd(f3,1,2) %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值

for n=2:1000000 %求等分数n

Rn=-(b-a)/12*((b-a)/n)^2*f2(x1) %计算余项

if abs(Rn)

break% 符合要求时结束

end

end

h=(b-a)/n %求h

Tn1=0

for k=1:n-1 %求连加和

xk=a+k*h

Tn1=Tn1+f(xk)

end

Tn=h/2*((f(a)+2*Tn1+f(b)))

z=exp(2)

R=Tn-z %求已知值与计算值的差

stem(xk,Tn1);

fprintf('用复化梯形算法计算的结果 Tn=')

disp(Tn)

fprintf('等分数 n=')

disp(n) %输出等分数

fprintf('已知值与计算值的误差 R=')

disp(R)

（2）复化Simpson

clc

clear

syms x%定义自变量x

f=inline('x^6/10-x^2+x','x') %定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可

f2=inline('36*x^2','x') %定义f(x)的四阶导数，输入程序1里求出的f2即可

f3='-(36*x^2)'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，一边求最大值e=5*10^(-8) %精度要求值

a=0 %积分下限

b=2 %积分上限

x1=fminbnd(f3,1,2) %求负的四阶导数的最小值点，也就是求四阶导数的最大值点对应的x值

for n=2:1000000 %求等分数n

Rn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1) %计算余项

if abs(Rn)

break% 符合要求时结束

end

end