数值分析上机答案 Microsoft Word 文档Word文件下载.docx
《数值分析上机答案 Microsoft Word 文档Word文件下载.docx》由会员分享，可在线阅读，更多相关《数值分析上机答案 Microsoft Word 文档Word文件下载.docx（46页珍藏版）》请在冰豆网上搜索。
fplot方法
y='
x^2*sin（x^2-x-2）'
;
fplot（y,[-22]）
（ii）参数方法
t=linspace（0,2*pi,100）;
x=2*cos（t）;
y=3*sin（t）;
（iii）
x=-3:
0.1:
3;
y=x;
[x,y]=meshgrid（x,y）;
z=x.^2+y.^2;
surf（x,y,z）
（iv）
y=-3:
13;
z=x.^4+3*x.^2+y.^2-2*x-2*y-2*x.^2.*y+6;
（v）
t=0:
x=sin（t）;
y=cos（t）;
z=cos（2*t）;
plot3（x,y,z）
（vi）
theta=linspace（0,2*pi,50）;
fai=linspace（0,pi/2,20）;
[theta,fai]=meshgrid（theta,fai）;
x=2*sin（fai）.*cos（theta）;
y=2*sin（fai）.*sin（theta）;
z=2*cos（fai）;
（vii）
x=linspace（0,pi,100）;
y1=sin（x）;
y2=sin（x）.*sin（10*x）;
y3=-sin（x）;
plot（x,y1,x,y2,x,y3）
实验5
程序exa_5
functionz=exa_5（x,y）
%thesizeofxandymustbesame.
z=zeros（size（x））;
[x1,y1]=find（x+y>
1）;
z（x1,y1）=0.5457*exp（-0.75*y（x1,y1）.^2-3.75*x（x1,y1）.^2-1.5*x（x1,y1））;
[x2,y2]=find（x+y<
=1&
x+y>
-1）;
z（x2,y2）=0.7575*exp（-y（x2,y2）.^2-6*x（x2,y2）.^2）;
[x3,y3]=find（x+y<
=-1）;
z（x3,y3）=0.5457*exp（-0.75*y（x3,y3）.^2-3.75*x（x3,y3）.^2+1.5*x（x3,y3））;
z=exa_5（x,y）;
实验6
查询trapz的功能和用法
helptrapz
查找trapz.m文件所在目录
whichtrapz
查看trapz.m程序结构
typetrapz
查看trapz.m文件所在目录还有哪些文件
cdc:
\MATLAB6p5\toolbox\matlab\datafun
dir
附录b上机实验题
实验1
symsthetafai
y=sin（fai）*cos（theta）-cos（fai）*sin（theta）;
simple（y）
ans=sin（fai-theta）
symsx;
f=x^4-5*x^3+5*x^2+5*x-6;
factor（f）
ans=（x-1）*（x-2）*（x-3）*（x+1）
symsa
A=[12;
2a]
在命令行中运行（求逆）
inv（A）
ans=
[a/（a-4）,-2/（a-4）]
[-2/（a-4）,1/（a-4）]
在命令行中运行（求特征值）
eig（A）
[1/2+1/2*a+1/2*（17-2*a+a^2）^（1/2）]
[1/2+1/2*a-1/2*（17-2*a+a^2）^（1/2）]
（1）
symskn
symsum（k^2,k,1,n）
ans=1/3*（n+1）^3-1/2*（n+1）^2+1/6*n+1/6
（2）
symsk
symsum（1/k^2,k,1,inf）
ans=1/6*pi^2
（3）
symsxn
symsum（1/（2*n+1）/（2*x+1）^（2*n+1）,n,0,inf）
simple（ans）（化简）
ans=1/2*log（（（（2*x+1）^2）^（1/2）+1）/（（（2*x+1）^2）^（1/2）-1））
symsxyz
f=sin（x^2*y*z）;
s=diff（diff（f,x,2）,y）;
s=subs（s,x,1）;
s=subs（s,y,1）;
s=subs（s,z,3）
s=88.2784
symsx
taylor（exp（x）,9）
taylor（log（1+x）,9）
taylor（sin（x）,9）
（4）
taylor（tan（x）,9）
实验7
symsy
int（exp（2*y）/（exp（y）+2））
diff（ans）%验证
simple（ans）%化简
symsxa
f=x^2/sqrt（a^2-x^2）;
int（f,x）
simple（diff（ans））%验证
symsabx
f=1/x/（sqrt（log（x）+a）+sqrt（log（x）+b））;
实验8见各章的上机答案
实验9
程序exb_9_11.m
s=1;
forn=2:
100
s=s*n;
s
tic;
exb_9_11;
toc
s=9.3326e+157
elapsed_time=0
程序exb_9_12.m
a=sym（1:
100）;
s=prod（a）;
exb_9_12;
s=93326……
elapsed_time=0.2500%可用double（s）转化为数值型
符号运算
程序
788657……%其375位数
数值计算的运算结果为Inf
第一章上机实验题
程序ex1_1.m
s=0;
forn=1:
1000
s=0.1+s;
s=s-100
ex1_1
s=-1.4069e-012
不稳定算法
程序ex1_2_1
a=zeros（1,21）;
a
（1）=1-1/exp
（1）;
21
a（n）=1-（n-1）*a（n-1）;
a
ex1_2_1
0.63210.36790.26420.20730.17090.14550.12680.11240.10090.0916
0.08390.07740.07180.06690.06270.05900.05550.0572-0.02951.5596
-30.1924
稳定算法
程序ex1_2_2
%ex1_2_2.m
a（21）=0.5*（1/exp
（1）+1）/21;
forn=20:
-1:
1
a（n）=（1-a（n+1））/n;
ex1_2_2
0.08390.07740.07180.06690.06270.05900.05570.05280.05010.0484
0.0326
a=[1.36110.75000.5250
0.75000.43600.3000
0.52500.30000.2136];
b=[2.6361;
1.4736;
1.0386];
a\b
1.1788
-0.0006
1.9658
程序ex1_3
functiony=ex1_4（p,x）
np=length（p）;
y=zeros（size（x））;
np
y=x.*y+p（n）;
p=[1030-26];
x=[1.11.21.3];
ex1_4（p,x）
9.403511.272313.7039
可用命令polyval（p,x）进行检验
第二章上机实验题
a=[41-1;
32-6;
1-53];
b=[9;
-2;
1];
det（a）
定解方程组
唯一解
a=[4-33;
b=[-1-21]'
a=[41;
32;
1-5];
b=[111]'
[rank（a）rank（[ab]）]
矛盾方程组
最小二乘解
a=[21-11;
121-1;
1121];
b=[123]'
null（a）
不定解方程组
无穷多解
det（a）%计算行列式
inv（a）%计算逆
[vd]=eig（a）%计算特征值及特征向量
[norm（a,1）norm（a,2）norm（a,inf）]%计算各种常用范数
[cond（a,1）cond（a,2）cond（a,inf）]%计算各种常用条件数
与
（1）类似
程序ex2_2_4.m
functiona=ex2_2_4（n）
a=zeros（n）;
a（1,[12]）=[56];
n-1
a（i,[i-1ii+1]）=[156];
a（n,[n-1n]）=[15];
命令
a=ex2_2_4（5）
a=ex2_2_4（50）
a=ex2_2_4（500）
[l,u,p]=lu（a）;
inv（p）*l*u
eig（a）
不是正定矩阵
a=[431;
33-5;
a=[5765;
71087;
68109;
57910];
r=chol（a）;
r'
*r
正定矩阵
不是对称矩阵
a=hilb（12）;
cond（a）
inv（a）
程序ex2_7.m
a=[11-1;
12-2;
-211];
b=[101]'
[l,u]=nalu（a）;
x=ex2_7l（l,b）;
x=ex2_7u（u,x）
程序ex2_7l.m
functionx=ex2_7l（l,b）
n=length（b）;
x=zeros（n,1）;
x
（1）=b
（1）/l（1,1）;
fork=2:
x（k）=（b（k）-l（k,1:
k-1）*x（1:
k-1））/l（k,k）;
程序ex2_7u.m
functionx=ex2_7u（u,b）
x（n）=b（n）/u（n,n）;
fork=（n-1）:
x（k）=（b（k）-u（k,k+1:
n）*x（k+1:
n））/u（k,k）;
实验8
a=[0.3e-1559.1431;
5.291-6.13-12;
11.2952;
1211];
b=[59.1746.7812]'
x=nagauss（a,b）
x=nagauss2（a,b）
1.2;
[a2,b2]=ex2_9_1（x）;
cond（a2）
1.5;
[a5,b5]=ex2_9_1（x）;
cond（a5）
1.8;
[a8,b8]=ex2_9_1（x）;
cond（a8）
程序ex2_9_1.m
function[a,b]=ex2_9_1（x）
n=length（x）-1;
x=x（:
a=ones（n+1,n+1）;
fork=1:
a（:
k+1）=x.^k;
b=sum（a'
）'
a5（6,6）=a5（6,6）+1e-4;
a5\b5
实验10
[a,b]=ex2_10_1（300）;
x=nagauss2（a,b,1）,toc
程序ex2_10_1.m
function[a,b]=ex2_10_1（n）
a（1,[12]）=[21];
a（k,k-1:
k+1）=[121];
a（n,[n-1,n]）=[12];
tic,x=ex2_10_2（a,b）;
程序ex2_10_2.m
functionx=ex2_10_2（a,b）
a（k,k）=a（k,k）-a（k,k-1）/a（k-1,k-1）*a（k-1,k）;
b（k）=b（k）-a（k,k-1）/a（k-1,k-1）*b（k-1）;
x（n）=b（n）/a（n,n）;
fork=n-1:
x（k）=（b（k）-a（k,k+1）*x（k+1））/a（k,k）;
x=a\b,toc
第三章上机实验题
p=[111];
x=roots（p）
polyval（p,x）
p=[30-402-1];
p=zeros（1,24）;
p（1,3）=-5;
p（1,7）=8;
p（1,8）=-6;
p（1,24）=5;
p=p（:
24:
x=roots（p）;
寻找有根区间或初始值
f=inline（'
x*log（sqrt（x^2-1）+x）-sqrt（x^2-1）-0.5*x'
fplot（f,[1,10]）;
gridon
有根区间为[1.52.5];
初始值为2
二分法
x=nabisect（f,1.5,2.5）
x=2.1155
Newton迭代法
df=inline（'
log（（x^2-1）^（1/2）+x）-0.5'
x=nanewton（f,df,2）
x^4-2^x'
fplot（f,[-22]）,gridon
[fzero（f,-0.8）,fzero（f,1.2）]
结果
ans=-0.86131.2396
[9*x
（1）^2+36*x
（2）^2+4*x（3）^2-36,x
（1）^2-2*x
（2）^2-20*x（3）,16*x
（1）-x
（1）^3-2*x
（2）^2-16*x（3）^2]'
'
x'
[x,y,h]=fsolve（f,[000]）
x=0.13420.9972-0.0985
clear
x1=sqrt（5）*cos（t）+2;
y1=sqrt（5）*sin（t）+3-2*x1;
x2=sqrt
（2）*cos（t）+3;
y2=6*sin（t）;
plot（x1,y1,x2,y2）
图形
显示结果
[（x
（1）-2）^2+（x
（2）-3+2*x
（1））^2-5,2*（x
（1）-3）^2+（x
（2）/3）^2-4]'
[x,y,h]=fsolve（f,[1.7-2.8]）
[x,y,h]=fsolve（f,[3.5-5.8]）
[x,y,h]=fsolve（f,[4-4]）
[x,y,h]=fsolve（f,[1.62]）
1.7362-2.6929
3.4829-5.6394
4.0287-4.1171
1.65811.8936
程序ex3_6.m
functionx=ex3_6（fname,x0,x1,eps,n）
ifnargin<
5n=50;
end
4eps=1e-6;
k=0;
whileabs（x0-x1）>
eps&
k<
=n
x=x1-feval（fname,x1）/（feval（fname,x1）-feval（fname,x0））*（x1-x0）;
x0=x1;
x1=x;
k=k+1;
x^3-2*x-5'
x=ex3_6（f,2.5,2）
x=2.0946
程序ex3_7.m
functionx=ex3_7（a,b,x0,eps,N）
5N=50;
3x0=zeros（n,1）;
d=diag（1./diag（a））;
lu=a-diag（diag（a））;
g=-d*lu;
f=d*b;
x=x0;
x0=x+2*eps;
whilenorm（x-x0）>
N
x0=x;
x=g*x0+f;
ifk==Nwarning（'
toomanytimes'
a=[521;
-132;
2-34];
b=[-1217-9]'
ex3_7（a,b,[000]'
0.5*1e-2）
ans=-4.0010
2.9965
2.0014
程序ex3_8.m
function[a,b]=ex3_8（n）
a（1,1:
3）=[3-1/2-1/4];
a（2,1:
4）=[-1/23-1/2-1/4];
a（n-1,n-3:
n）=[-1/4-1/23-1/2];
a（n,n-2:
n）=[-1/4-1/23];
n-2
a（i,i-2:
i+2）=[-1/4-1/23-1/2-1/4];
[a,b]=ex3_8（100）;
na
升级会员 