计算数值方法上机作业

作业一
4-11（1）用Romberg 方法求下列积分，要求误差不超过a
,8
.00
dx x a=10-2
计算软件MATLAB
程序代码：
第一步：建立M 文件
function [A,err]=romberg(f,a1,a2,e)
%A 为计算结果矩阵；err 为计算误差；
% romberg 参数依次为要积分的函数，下界，上界，精度；
h=a2-a1;
A(1,1)=(h/2)*(feval(f,a1)+feval(f,a2));
A(2,1)=A(1,1)/2+(h/2)*feval(f,a1+(h/2));
i=2;j=1;n=2;
while 1
m=1;
h=h/2;
while j<i
j=j+1;
A(i,j)=(4^m/((4^m)-1))*A(i,j-1)-A(i-1,j-1)/((4^m)-1);
m=m+1;
end
if abs(A(j,j)-A(j-1,j-1))<e
err=(abs(A(j,j)-A(j-1,j-1)))
break;
end
j=1;
i=i+1;
A(i,j)=0;
for k=1:n
A(i,j)=A(i,j)+feval(f,a1+(2*k-1)*(h/2));
end
A(i,j)=A(i,j)*(h/2);
A(i,j)=A(i-1,j)/2+A(i,j);
n=2*n;
end
首先在窗口中输入format long命令，令计算结果显示15位有效数字
>> f=inline('sqrt(x)')
f =
Inline function:
f(x) = sqrt(x)
>> romberg(f,0,0.8,1e-2)
err =
0.00418661034609
ans =
0.35777087639997 0 0 0
0.43186765101345 0.45656657588462 0 0
0.46029587845502 0.46977195426887 0.47065231282782 0
0.47091965235177 0.47446091031736 0.47477350738725 0.47483892317391
用Romberg 方法求得积分结果为0.47483892317391
作业二
计算软件：MATLAB
（1）Jacobi方法
程序代码：
第一步：建立M文件
function [x,k,flag,err]=Jacobi(A,b,delta,max1)
% A为方程组的系数矩阵；
% b为方程组的右端项
% delta为精度要求，默认为le-5;
% max1为最大迭代次数，默认为100；
% x为方程组的解；
% k为迭代次数；
% flag为指标变量flag=‘OK!’表示迭代收敛到指标要求；% flag=‘fail’表示迭代失败;
if nargin<4 max1=100;end
if nargin<3 delta=1e-5;end
n=length(A);k=0;
x=zeros(n,1);y=zeros(n,1);flag='OK!';
while 1
for i=1:n
y(i)=b(i);
for j=1:n
if j~=i
y(i)=y(i)-A(i,j)*x(j);
end
end
if abs(A(i,j))<1e-10|k==max1
flag='Fail';return;
end
y(i)=y(i)/A(i,i);
end
if norm(y-x,inf)<delta
return
end
x=y;k=k+1;
end
首先在窗口中输入format long命令，令计算结果显示15位有效数字然后在命令窗口输入系数矩阵和右端项：
>>A =[5 2 1;-1 4 2;2 -3 10];b=[-12 20 3]
回车得到：
A =
5 2 1
-1 4 2
2 -
3 10
b =
-12 20 3
再输入：
>> [x,k,flag]=Jacobi(A,b)
x =
-3.99998953291871
3.00000094812013
1.99999144958428
k =
19
flag =
OK!
这说明经过19次迭代得到满足精度要求的近似解
x=[-3.999 989 532 918 71, 3.000 000 948 120 13, 1.999 991 449 584 28]T
（2）Gauss-Seidel方法
程序代码：
第一步：建立M文件
function [x,k,flag,err]= Gau_Seid (A,b,delta,max1)
% A为方程组的系数矩阵；
% b为方程组的右端项
% delta为精度要求，默认为le-5;
% max1为最大迭代次数，默认为100；
% x为方程组的解；
% k为迭代次数；
% flag为指标变量flag=‘OK!’表示迭代收敛到指标要求；% flag=‘fail’表示迭代失败
if nargin<4 max1=100;end
if nargin<3 delta=1e-5;end
n=length(A);k=0;
x=zeros(n,1);y=zeros(n,1);flag='OK!';
while 1
for i=1:n
y(i)=b(i);
for j=1:n
if j~=i
y(i)=y(i)-A(i,j)*x(j);
end
end
if abs(A(i,j))<1e-10|k==max1
flag='Fail';return;
end
y(i)=y(i)/A(i,i);
end
if norm(y-x,inf)<delta
return
end
x=y;k=k+1;
end
然后在命令窗口输入系数矩阵和右端项：
>>A =[5 2 1;-1 4 2;2 -3 10];b=[-12 20 3]
回车得到：
A =
5 2 1
-1 4 2
2 -
3 10
b =
-12 20 3
再输入：
>> [x,k,flag]= Gau_Seid (A,b)
回车得到：
x =
-4.00000011851181
3.00000049280490
2.00000017154383
k =
9
flag =
OK!
这说明经过9次迭代得到满足精度要求的近似解
x=[-4.000 000 118 511 81, 3.000 000 492 804 90, 2.000 000 171 543 83]T。