数值分析试题.docx
- 文档编号:9734471
- 上传时间:2023-02-06
- 格式:DOCX
- 页数:14
- 大小:58.91KB
数值分析试题.docx
《数值分析试题.docx》由会员分享,可在线阅读,更多相关《数值分析试题.docx(14页珍藏版)》请在冰豆网上搜索。
数值分析试题
1.设有某实验数据如下:
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
5.1234
5.3053
5.5684
5.9378
6.4270
7.0798
7.9493
9.0253
10.3627
(1)在MATLAB中作图观察离散点的结构,用最小二乘法拟合一个合适的多项式函数;
(2)在MATLAB中作出拟合曲线图.
解:
(1)在MATLAB中作图观察离散点的结构,用最小二乘法拟合一个合适的多项式函数。
具体程序如下:
>>x=[0.10.20.30.40.50.60.70.80.9];
>>y=[5.12345.30535.56845.93786.42707.07987.94939.025310.3627];
>>plot(x,y,'r+'),legend('x,y'),xlabel('x'),ylabel('y'),title('(x,y)的散点结构图')
从图中各点可以看到,(x,y)散点结构图的变化趋势与二次多项式很接近,故选取曲线函数为二次多项式函数,拟合多项式的次数为2。
用最小二乘法拟合一个合适的多项式函数,具体程序如下:
>>s=polyfit(x,y,2);
>>P=poly2str(s,'t')
P=
8.2191t^2-1.8822t+5.3139
则曲线函数为发f(x)=8.2192t2-1.8822t+5.3139
(2)在MATLAB中作出拟合曲线图,具体程序如下:
>>x=linspace(0.1,0.9,9);
>>y=8.2191.*x.^2-1.8822.*x+5.3139;
>>plot(x,y),legend('x,y'),xlabel('x'),ylabel('y'),title('拟合曲线y=f(x)')
2.在MATLAB中用复合Simpson公式编程计算下列积分,并使其结果至少有6位有效数字.
解:
(1)所安装的MATLAB中不含有复合辛普森程序,用M文件编辑器编写被积函数的M文件,并以jfSimpson.m命名保存至MATLAB搜索路径下。
程序如下:
function[I,step]=jfSimpson(f,a,b,type,eps)
%type=1辛普森公式
%type=2复合辛普森公式
if(type==2&&nargin==8)
eps=1.0e-8;%缺省精度为0.0001
end
I=0;
switchtype
case1,
I=((b-a)/6)*(subs(sym(f),findsym(sym(f)),a)+...
4*subs(sym(f),findsym(sym(f)),(a+b)/2)+...
subs(sym(f),findsym(sym(f)),b));
step=1;
case2,
n=2;
h=(b-a)/2;
I1=0;
I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;
whileabs(I2-I1)>eps
n=n+1;
h=(b-a)/n;
I1=I2;
I2=0;
fori=0:
n-1
x=a+h*i;
x1=x+h;
I2=I2+(h/6)*(subs(sym(f),findsym(sym(f)),x)+...
4*subs(sym(f),findsym(sym(f)),(x+x1)/2)+...
subs(sym(f),findsym(sym(f)),x1));
end
end
I=I2;
step=n;
end
(2)在命令窗口中,调用jfSimpson函数进行求解。
命令窗口显示如下:
>>[I,s]=jfSimpson(inline('exp(x.^2)'),0,2,2,1.0e-8)
I=
16.4526
s=
102
说明:
求解积分函数为16.4526。
3.在MATLAB中试用经典四阶Runge-Kutta法编程求解初值问题
解:
(1)用M文件编辑器将RK4.m、RK_fun.m、RK_main.m的主程序命名并保存至MATLAB搜索路径下。
主程序如下:
RK4.m
function[t,y]=RK4(func,t0,tt,y0,N,varagin)
%Rk方法计算一阶级微分方程组,
%由微分方程知识可以知道,对于高阶微分方程可以化为一阶微分方程组。
%初始时刻为t0,结束时为tt,初始时刻函数值为y0
%N为步数
ifnargin<4
N=100
end
%如果没有输入N的话,那么默认计算步长为(tt-t0)/100
y(1,:
)=y0(:
)';
h=(tt-t0)/N;%步长
t=t0+[0:
N]'*h;
fork=1:
N
f1=h*feval(func,t(k),y(k,:
));
f1=f1(:
)';
%RK中的k1
f2=h*feval(func,t(k)+h/2,y(k,:
)+f1/2);
f2=f2(:
)';
%RK中的k2
f3=h*feval(func,t(k)+h/2,y(k,:
)+f2/2);
f3=f3(:
)';
%RK中的k3
f4=h*feval(func,t(k)+h,y(k,:
)+f3);
f4=f4(:
)';
%RK中的k4
y(k+1,:
)=y(k,:
)+(f1+2*(f2+f3)+f4)/6;
%所求结果
end
RK_fun.m
%RK4方法的测试函数
functionf=RK_fun(t,y)
f=-3*y^2+2*t^2+3*t;
RK_main.m
%rk方法的主函数
x0=0;
xt=1;
y0=1;
%符号计算给出分析解
y=dsolve('Dy=-3*y^2+2*t^2+3*t','y(0)=1','t')
%RK4方法给出计算结果
[x,yrk]=RK4('RK_fun',x0,xt,y0,10);
%对应于真解的离散数据
yreal=subs(y,x);
tol=yreal-yrk;%RK4方法与分析解的误差
plot(x,yreal,x,yrk,'+',x,tol,'*')
legend('分析解','RK4计算结果','两者误差')
yrk;
[x,yrk,tol]
%yrk为所要计算的结果
(2)运行结果如下:
>>RK_main
y=
<1x1sym>
ans=
<1x1sym><1x1sym><1x1sym>
<1x1sym><1x1sym><1x1sym>
<1x1sym><1x1sym><1x1sym>
<1x1sym><1x1sym><1x1sym>
<1x1sym><1x1sym><1x1sym>
<1x1sym><1x1sym><1x1sym>
<1x1sym><1x1sym><1x1sym>
<1x1sym><1x1sym><1x1sym>
<1x1sym><1x1sym><1x1sym>
<1x1sym><1x1sym><1x1sym>
<1x1sym><1x1sym><1x1sym>
4.在MATLAB中利用Newton法编程计算方程
的一个根.
解:
(1)所安装的MATLAB中不含有newton.m程序,用M文件编辑器编写被积函数的M文件,并以newtonqxf.m命名保存至MATLAB搜索路径下。
程序如下:
functionroot=newtonqxf(f,a,b,eps)
%牛顿法求函f数在区间[a,b]上的一个零点
%f为函数名
%a为区间左端点
%b为区间右端点
%eps为根的精度
%root为求出的函数零点
if(nargin==3)
eps=1.0e-4;
end
f1=subs(sym(f),findsym(sym(f)),a);
f2=subs(sym(f),findsym(sym(f)),b);
if(f1==0)
root=a;
end
if(f2==0)
root=b;
end
if(f1*f2>0)
disp('两端点函数值乘积大于0!
');
return;
else
tol=1;
fun=diff(sym(f));
fa=subs(sym(f),findsym(sym(f)),a);
fb=subs(sym(f),findsym(sym(f)),b);
dfa=subs(sym(fun),findsym(sym(fun)),a);
dfb=subs(sym(fun),findsym(sym(fun)),b);
if(dfa>dfb)
root=a-fa/dfa;
else
root=b-fb/dfb;
end
while(tol>eps)
r1=root;
fx=subs(sym(f),findsym(sym(f)),r1);
dfx=subs(sym(fun),findsym(sym(fun)),r1);
root=r1-fx/dfx;
tol=abs(root-r1);
end
end
(2)在命令窗口中,调用newtonqxf函数进行求解。
命令窗口显示如下:
>>root=newtonqxf('x-log(x)-2',3,1.0e-9)
root=
3.1462
5.已知线性方程组
(1)利用MATLAB说明解该方程组的Gauss-Seidel迭代法是否收敛;
(2)在MATLAB中利用Gauss-Seidel迭代法求解该方程组,并给出迭代次数.
解:
通过M文件编辑器加入gauseidel.m函数程序,如下:
function[x,n]=gauseidel(A,b,x0,eps,M)
%求解线性方程组的迭代法,其中,
%A为方程组的系数矩阵;
%b为方程组的右端项;
%x0为迭代初始化向量
%eps为精度要求,缺省值为1e-5;
%M为最大迭代次数,缺省值100;
%x为方程组的解;
%n为迭代次数;
ifnargin==3
eps=1.0e-6;
M=200;
elseifnargin==4
M=200;
elseifnargin<3
error
return;
end
D=diag(diag(A));%求A的对角矩阵
L=-tril(A,-1);%求A的下三角阵
U=-triu(A,1);%求A的上三角阵
G=(D-L)\U;
f=(D-L)\b;
x=G*x0+f;
n=1;%迭代次数
whilenorm(x-x0)>=eps
x0=x;
x=G*x0+f;
n=n+1;
if(n>=M)
disp('Warning:
迭代次数太多,可能不收敛!
');
return;
end
end
(1)在命令窗口中,输入程序如下:
>>A=[5,2,1;-1,4,2;2,-5,10]
A=
521
-142
2-510
>>D=diag(diag(A))
D=
500
040
0010
>>L=-tril(A,1)
L=
-5-20
1-4-2
-25-10
>>U=-triu(A,1)
U=
0-2-1
00-2
000
>>G=(D-L)\U
G=
0-0.1945-0.0515
0-0.0275-0.2426
00.0126-0.0555
>>p=abs(eigs(G,1))
p=
0.0677
因p=0.0677<1,该矩阵方程组Gauss-Seidel迭代收敛。
(2)命令窗口中输入程序如下:
>>A=[5,2,1;-1,4,2;2,-5,10];
>>b=[-12,10,1]';
>>x0=zeros(3,1);
>>[x,n]=gauseidel(A,b,x0)
x=
-3.0909
1.0945
1.2655
n=
11
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 试题