数值分析王能超matlab上机作业.docx
《数值分析王能超matlab上机作业.docx》由会员分享，可在线阅读，更多相关《数值分析王能超matlab上机作业.docx（16页珍藏版）》请在冰点文库上搜索。
数值分析王能超matlab上机作业
引论习题
题目：
用二分法求方程x^3-x-1=0在[1,2]内的近似根，要求误差不超过10-3
程序：
function[root,n]=dichotomy（a,b,eps）
%dichotomy：
二分法求根函数
%f:
带求解方程；[a,b]求解区间；n为二分次数
ifnargin<2disp（'输入变量不够'）;return;
elseifnargin>3disp（'输入变量过多'）;return;
elsenargin==2
eps=1.0e-3;
end
end
ifa>bdisp（'区间输入错误'）
return
end
c=（a+b）/2;
iff（c）==0
y=c
n=1
return
end
n=0;
whileabs（b-a）>eps
iff（a）*f（c）>0
a=c;
elseb=c;
end
c=（a+b）/2;
n=n+1;
end
root=c;
end
运算结果：
定义函数f如下：
function[yy]=f（x）
%定义了测试函数
%x为输入变量
yy=x^3-x-1;
end
在commondwindow中运行结果如下：
>>[root,n]=dichotomy（1,2,1.0e-3）
root=
1.32471
n=
10
习题一
题目一给出概率积分
的数据表
i
0
1
2
3
xi
0.46
0.47
0.48
0.49
yi
0.4846555
0.4937452
0.5027498
0.5116683
用二次插值计算，试问：
（1）当x=0.472时该积分值等于多少？
（2）当x为何值时积分值等于0.5？
程序：
function[y,t]=interpolation（X,Y,x）
%拉格朗日插值函数
%X为自变量数组，Y为函数值数组，x为插值点，t为插值基函数
n=length（X）;
y=0;t=zeros（1,n）;
fori=1:
n
t（i）=1;
forj=1:
n
ifi==jcontinue
end
t（i）=t（i）*（x-X（j））/（X（i）-X（j））;
end
y=y+t（i）*Y（i）;
end
end
运算结果：
在commondwindow中运行结果如下：
（1）当x=0.472时该积分值等于多少？
>>clear
>>X=[0.460.470.480.49];
>>Y=[0.48465550.49374520.50274980.5116683];
>>x=0.472;
>>[y,t]=interpolation（X,Y,x）
y=
0.4956
t=
-0.04800.86400.2160-0.0320
（2）当x为何值时积分值等于0.5？
>>clear
>>X=[0.460.470.480.49];
>>Y=[0.48465550.49374520.50274980.5116683];
>>y=0.5;
>>[x,t]=interpolation（Y,X,y）
x=
0.4769
t=
-0.04520.33560.7707-0.0611
题目二构造适合下列数据表的三次插值样条函数S（x）:
x
-1
0
1
3
y
-1
1
3
5
y'
6
1
程序：
function[y,m]=spline（X,Y,x,m1,mn）
%样条插值函数
%X，Y为输入的自变量、函数值数组；
%x为插值点（数组），y为插值结果（数组），m为各插值节点的一阶导数值
%m1为X
（1）的一阶导数;mn为X（n）的一阶导数
n=length（X）;
alpha=zeros（n-2,1）;
beta=zeros（n-2,1）;
fori=1:
（n-2）
alpha（i）=（X（i+1）-X（i））/（X（i+1）-X（i）+X（i+2）-X（i+1））;
beta（i）=3*（（1-alpha（i））*（Y（i+1）-Y（i））/（X（i+1）-X（i））+alpha（i）*（Y（i+2）-Y（i+1））/（X（i+2）-X（i+1）））;
end
m=zeros（n,1）;
m
（1）=m1;m（n）=mn;
A=zeros（（n-2）,（n-2））;
B=zeros（（n-2）,1）;
forj=1:
（n-2）
A（j,j）=2;
end
fork=2:
（n-2）
A（k,（k-1））=1-alpha（k）;
end
forl=1:
（n-3）
A（l,（l+1））=alpha（l）;
end
B
（1）=beta
（1）-（1-alpha
（1））*m1;
B（n-2）=beta（n-2）-alpha（n-2）*mn;
forp=2:
（n-3）
B（p）=beta（p）;
end
m（2:
（n-1））=A\B;
s=length（x）;
y=zeros（1,s）;
fort=1:
s
r=1;
forq=1:
n
ifx（t）>=X（q）&&x（t）<=X（q+1）
break;
end
r=r+1;
end
xx=（x（t）-X（r））/（X（r+1）-X（r））;
y（t）=（xx-1）^2*（2*xx+1）*Y（r）+xx^2*（-2*xx+3）*Y（r+1）+（X（r+1）-X（r））*xx*（xx-1）^2*m（r）+（X（r+1）-X（r））*xx^2*（xx-1）*m（r+1）;
end
end
运算结果：
在commondwindow中运行结果如下：
>>clear
>>X=[-1013];
>>Y=[-1135];
>>m1=6;mn=1;
>>x=-1:
0.1:
3;
>>[y,m]=spline（X,Y,x,m1,mn）;
plot（x,y）
求得结果如下（整理之后）：
x
-1
-0.9
-0.8
-0.7
-0.6
-0.5
-0.4
-0.3
-0.2
-0.1
y
-1
-0.47
-0.056
0.251
0.472
0.625
0.728
0.799
0.856
0.917
x
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
y
1
1.119
1.272
1.453
1.656
1.875
2.104
2.337
2.568
2.791
x
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
y
3
3.19
3.362
3.517
3.656
3.781
3.894
3.996
4.088
4.172
x
2
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
y
4.25
4.323
4.392
4.459
4.526
4.594
4.664
4.738
4.818
4.905
x
3
y
5
其中红色数据表示插值节点
将以上结果绘制成图，如下所示:
习题二
题目编写一通用型复化辛普生公式，能够对任意长度的等间距离散数据进行积分运算。
要求：
命令格式：
y=simpson（x,h）;
y：
积分结果
x：
输入数组
h：
间距
程序：
function[y]=simpson（x,h）
%一通用型复化辛普森求积公式
%y:
积分结果；x:
输入数组；h:
间距
ifnargin<2disp（'输入参数不够，请准确输入积分数据和步长'）;return;
end
n=length（x）;
%X=zeros（n+1）;
ifn==1y=x;return;
elseifmod（n,2）==1
%X=zeros（n+1）;
A=zeros（（n+1）/2）;
B=zeros（（n-1）/2）;
A=x（1:
2:
（end-2））;
B=x（2:
2:
（end-1））;
y=2*h/6*（x（end）-x
（1）+4*sum（B）+2*sum（A））;
%X（1:
n）=x;
elsemod（n,2）==1
X=zeros（n-1）;
X=x（1:
（end-1））;
m=length（X）;
A=zeros（（m+1）/2）;
B=zeros（（m-1）/2）;
A=X（1:
2:
（end-2））;
B=X（2:
2:
（end-1））;
y=2*h/6*（X（end）-X
（1）+4*sum（B）+2*sum（A））+（x（n）+x（n-1））/2*h;
end
end
end
运算结果：
以计算
为例
在commondwindow中运行结果如下：
clear
>>t=0:
（1/8）:
1;
>>x=sin（t）./t;
>>h=1/8;
>>[y]=simpson（x,h）
y=
0.9461
以上积分结果与教材P65所给结果完全一致
习题三
题目分别用显式和隐式的二阶亚当姆斯方法求解初值问题y'=1-y，y（0）=0，令y（0.2）=0.181，取h=0.2计算y（1.0）。
程序：
1、显式亚当姆斯方法：
function[X,Y]=Adams_shown（x0,y0,y1,h,N）
%显式二阶亚当姆斯求解微分方程
%x0，y0,y1为输入初值，h为步长，N为计算点数（包括除x0,y0在内的输入初值）
%X,Y分别为输出数组,f为y'=f（x,y）
X=zeros（N+1,1）;Y=zeros（N+1,1）;%为X,Y预置空间
X=linspace（x0,x0+N*h,N+1）;X=X（:
）;
Y
（1）=y0;Y
（2）=y1;
forn=2:
N
yy0=f（x0,y0）;yy1=f（x0+h,y1）;
y2=y1+h/2*（3*yy1-yy0）;
yy2=f（x0+2*h,y2）;
Y（n+1）=y2;
x0=x0+h;
y0=y1;
y1=y2;
end
end
2、隐式亚当姆斯方法：
function[X,Y]=Adams_hidden（x0,y0,h,N）
%隐式二阶亚当姆斯求解微分方程
%x0，y0为输入初值，h为步长，N为计算点数（不包括输入初值）
%X,Y分别为输出数组,f为y'=f（x,y）
X=zeros（N+1,1）;Y=zeros（N+1,1）;%为X,Y预置空间
X=linspace（x0,x0+N*h,N+1）;X=X（:
）;
Y
（1）=y0;
forn=1:
N
yy0=f（x0,y0）;
%y1=solve（'y0+h/2*（f（x0+h,y1）+yy0）-y1','y1'）;
y1=（y0+h/2*（1+yy0））/（1+h/2）;
yy1=f（x0+h,y1）;
Y（n+1）=y1;
x0=x0+h;
y0=y1;
end
end
3、二阶亚当姆斯预报—校正方法：
function[X,Y]=Adams_shown_hidden（x0,y0,y1,h,N）
%二阶亚当姆斯预报-校正求解微分方程
%x0，y0,y1为输入初值，h为步长，N为计算点数（包括除x0,y0在内的输入初值）
%X,Y分别为输出数组,f为y'=f（x,y）
X=zeros（N+1,1）;Y=zeros（N+1,1）;%为X,Y预置空间
X=linspace（x0,x0+N*h,N+1）;X=X（:
）;
Y
（1）=y0;Y
（2）=y1;
forn=2:
N
yy0=f（x0,y0）;yy1=f（x0+h,y1）;
y2=y1+h/2*（3*yy1-yy0）;
yy2=f（x0+2*h,y2）;
y2=y1+h/2*（yy2+yy1）;
yy2=f（x0+2*h,y2）;
Y（n+1）=y2;
x0=x0+h;
y0=y1;
y1=y2;
end
end
运算结果：
定义函数y'=1-y如下：
function[yy]=f（x,y）
%定义了测试函数
%x，y为输入变量
yy=1-y;
end
在commondwindow中运行结果如下：
>>clear
>>x0=0;y0=0;y1=0.181;h=0.2;N=5;
>>[X1,Y1]=Adams_shown（x0,y0,y1,h,N）;
>>[X2,Y2]=Adams_hidden（x0,y0,h,N）;
>>[X3,Y3]=Adams_shown_hidden（x0,y0,y1,h,N）;
>>plot（X1,Y1,'r',X2,Y2,'b',X3,Y3,'g'）;
>>legend（'显式亚当姆斯','隐式亚当姆斯',,'二阶亚当姆斯校正-预报方法'）
最后所得结果如下表所示：
X
0
0.2
0.4
0.6
0.8
1
显式
0
0.181
0.3267
0.4468
0.5454
0.6265
隐式
0
0.1818
0.3306
0.4523
0.5519
0.6334
校正-预报
0
0.181
0.3302
0.4523
0.5521
0.6337
根据上述数据绘制成图，如下所示：
习题四
题目用牛顿法求下列方程的根，要求计算结果小数点后有4位有效数字
（1）
（2）
程序：
functionroot=NewtonRoot（f,x0,eps）
%牛顿法求方程的根
%x0为迭代初值
%f为方程表达式：
f（x）=0
%eps为迭代精度
%root为满足精度要求的近似根
if（nargin==2）
eps=1.0e-4;
end
tol=1;
fx=sym（f）;
dfx=diff（sym（f））;%求导数
while（tol>eps）
fx0=subs（fx,findsym（fx）,x0）;
dfx0=subs（dfx,findsym（dfx）,x0）;%求该点的导数值
root=x0-fx0/dfx0;%迭代的核心公式
tol=abs（root-x0）;
x0=root;
end
formatlongg;
root=eval（root）;
end
运算结果：
在commondwindow中运行结果如下：
（1）
>>clear
>>root=NewtonRoot（'x^3-3*x-1',2,1.0e-4）
root=
1.87938524483667
（2）
>>root=NewtonRoot（'x^2-3*x-exp（x）+2',1,1.0e-4）
root=
0.257530285426349
升级会员 