计算机.docx
《计算机.docx》由会员分享，可在线阅读，更多相关《计算机.docx（15页珍藏版）》请在冰豆网上搜索。
计算机
《MATLAB程序设计实践》课程考核
班级：
材料0808
学号：
0604080805
姓名：
韩勇军
1、编程实现以下科学计算算法，并举一例应用之。
（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）
“切比雪夫逼近”
算法说明：
当一个连续函数定义在区间[-1,1]上时，它可以展开成切比雪夫级数。
即：
其中为次切比雪夫多项式，具体表达可通过递推得出：
，
它们之间满足如下的正交关系：
在实际应用中，可根据所需的精度来截取有限的项数，切比雪夫级数中的系数由下式决定：
在MATLAB中编程实现的切比雪夫逼近法函数为：
Chebyshev。
功能：
用切比雪夫多项式逼近已知函数。
调用格式：
其中，y为已知函数；
k为逼近已知函数所需项数；
f是求得的切比雪夫逼近多项式在x0处的逼近值。
程序源代码（m文件）：
functionf=Chebyshev（y,k,x0）
%用切比雪夫多项式逼近已知函数
%已知函数：
y
%逼近已知函数所需项数：
k
%逼近点的x坐标：
x0
%求得的切比雪夫逼近多项式或在x0处的逼近值：
f
symst;
T（1:
k+1）=t;
T
（1）=sym（'1'）;
T
（2）=t;
c（1:
k+1）=sym（'0'）;
c
（1）=int（subs（y,findsym（sym（y））,sym（'t'））*T
（1）/sqrt（1-t^2）,t,-1,1）/pi;
c
（2）=2*int（subs（y,findsym（sym（y））,sym（'t'））*T
（2）/sqrt（1-t^2）,t,-1,1）/pi;
f=c
（1）+c
（2）*t;
fori=3:
k+1
T（i）=2*t*T（i-1）-T（i-2）;
c（i）=2*int（subs（y,findsym（sym（y））,sym（'t'））*T（i）/sqrt（1-t^2）,t,-1,1）/2;
f=f+c（i）*T（i）;
f=vpa（f,6）;
if（i==k+1）
if（nargin==3）
f=subs（f,'t',x0）;
else
f=vpa（f,6）;
end
end
end
应用实例：
切比雪夫应用实例。
用切比雪夫公式（取6项）逼近函数，并求当x=0.5时的函数值。
解：
利用程序求解方程，在MATLAB命令窗口中输入：
>>Chebyshev（'1/（2-x）',6）%调用创建的函数euler，输出切比雪夫多项式的6个项
再在MATLAB命令窗口中输入：
>>Chebyshev（'1/（2-x）',6，0.5）%调用创建的函数euler，输出当x=0.5时的函数值
输出结果：
流程图：
2、编程解决以下科学计算问题。
实验2：
用二分法（程序3.1）和Newton迭代法（程序3.2）求下列方程的正根：
（1）二分法：
erfen.m:
functiony=erfen（fun,a,b,esp）
ifnargin<4esp=1e-4;end
iffeval（fun,a）*feval（fun,b）<0
n=1;
c=（a+b）/2;
whilec>esp
iffeval（fun,a）*feval（fun,c）<0
b=c;c=（a+b）/2;
elseiffeval（fun,c）*feval（fun,b）<0
a=c;c=（a+b）/2;
elsey=c;esp=10000;
end
n=n+1;
end
y=c;
elseiffeval（fun,a）==0
y=a
elseiffeval（fun,b）==0
y=b;
elsedisp（'these,maynotbearootintheintercal'）;
end
n
函数文件：
fc.m
functiony=fc（x）
y=x*log（sqrt（x.^2-1）+x）-sqrt（x.^2-1）-0.5*x;
end
在MATLAB命令框中输入：
>>clear
>>erfen（'fc',1,10,0.05）
输出结果：
（2）牛顿迭代法：
newton.m
functiony=newton（x0）
x1=x0-fc（x0）/df（x0）;
n=1;
while（abs（x1-x0）>=1.0e-6）&（n<=100000000）
x0=x1;
x1=x0-fc（x0）/df（x0）;n=n+1;
end
x1
n
函数文件：
fc.m
functiony=fc（x）
y=x*log（sqrt（x.^2-1）+x）-sqrt（x.^2-1）-0.5*x;
end
导数函数文件：
df.m
functiony=df（x）
y=log（x+（x^2-1）^（1/2））-x/（x^2-1）^（1/2）+（x*（x/（x^2-1）^（1/2）+1））/（x+（x^2-1）^（1/2））-1/2;
end
在MATLAB命令框中输入：
>>clear
>>clc
>>newton（3）
输出结果：
对比二分法和牛顿法可知:
newton法n值较小，收敛得要快一些。
实验3：
已知函数f（x）=x4-2x在（-2,2）内有两个根，作图求取初值。
用MATAB命令求这两个根。
源代码：
erfen2.m:
functiony=erfen2（fun,a,b,n）
formatlong
ifnargin<4
n=10000;
end
if（bdisp（'初始值错误！
'）
y=NaN;
return;
end;
if（feval（fun,a）*feval（fun,b）<0）
c=a/2+b/2;
while（n>0）
iffeval（fun,a）*feval（fun,c）<0
b=c;c=a/2+b/2;
elseiffeval（fun,b）*feval（fun,c）<0
a=c;c=a/2+b/2;
else
y=c;n=0;
end
n=n-1;
end
y=c;
elseiffeval（fun,a）==0
y=a;
elseiffeval（fun,b）==0;
y=b;
else
disp（'取值范围内可能不存在解！
'）
end
函数文件：
f2.m
functiony=f2（x）
y=x^4-2^x;
end
在MATLAB命令框中输入：
>>fplot（'f2',[-2,2]）
>>
得到函数图形：
从图形上可以看出f（-2）,f
（2）的值均大于0，而f（0）的值为小于。
下面对其进行验证：
在命令框中分别输入：
>>y=f2
（2）
>>
>>y=f2（-2）
>>
>>y=f2（0）
>>
得出的结果为：
该结果验证了上述猜想：
故易知函数的两个解分别位于区间[-2,0]，[0,2]中。
下面为具体的求解过程并分析误差大小：
在MATLAB命令框中输入：
>>x=erfen2（'f2',-2,0）
x=
-0.861345332309651
>>tol=f2（x）
tol=
2.220446049250313e-016
>>x=erfen2（'f2',0,2）
x=
1.239627729522762
>>tol=f2（x）
tol=
8.881784197001252e-016
输出结果：
在MATLAB中实现切比雪夫逼近的代码如下：
M文件:
程序：
r=Secant（'log（x）+sqrt（x）-2',1,4）
运行结果：
r=
1.8773
1、编程解决以下科学计算问题
计算积分
1
程序：
>>symsx
>>int（exp（-x.^2）./（1+x.^2）,-inf,+inf）
ans=
pi*erfc
（1）*exp
（1）
②
M文件：
functiony=fun（x）
y=tan（x）./x.^（0.7）;
程序：
quad（'fun',0,1）
结果：
ans=
0.9063
③
M文件：
functiony=fun1（x）
y=exp（x）./（1+x.^2）.^0.5;
程序：
quad（'fun1',0,1）
结果：
ans=
1.4690
④
程序及结果：
symsxy
>>int（1+x+y.^2,y,-sqrt（1-（x-1）.^2）,sqrt（1-（x-1）.^2））
ans=
（2*（1-（x-1）^2）^（1/2）*（3*x-（x-1）^2+4））/3
>>int（（2*（1-（x-1）^2）^（1/2）*（3*x-（x-1）^2+4））/3,x,-1,1）
ans=
pi/2+（5*asin
（2））/4+（3*3^（1/2）*i）/2-2/3+log（2-3^（1/2））*i
编制一个变步长梯行程序，解例题5.14比较与Romberg积分法的计算量
例5.14用Romberg积分公式计算积分I＝时，要求结果误差为0.5×10。
变步长梯行法：
M文件：
function[I,step]=BBC（f,a,b,eps）
if（nargin==3）
eps=0.5e-5;
end;
n=1;
h=（b-a）;
I2=（subs（sym（f）,findsym（sym（f））,a）+subs（sym（f）,findsym（sym（f））,b））*h/2;
Tol=1;
whileTol>eps
n=n+1;
h=h./2;
I1=I2;
I2=0;
fori=0:
（2^n-1）
x=a+h.*i;
x1=x+h;
I2=I2+（h./2）.*（subs（sym（f）,findsym（sym（f））,x）+subs（sym（f）,findsym（sym（f））,x1））;
end
Tol=abs（I2-I1）;
end
I=I2;
step=n;
Romberg积分法：
M文件：
function[I,step]=Roberg（f,a,b,eps）
if（nargin==3）
eps=0.5e-5
end;
symsf
M=1;
tol=10;
k=0;
T=zeros（1,1）;
h=b-a;
T（1,1）=（h/2）.*（subs（f,findsym（f）,a）+subs（f,findsym（f）,b））;
whiletol>eps
k=k+1;
h=h/2;
Q=0;
fori=1:
M
x=a+h*（2*i-1）;
Q=Q+subs（f,findsym（f）,x）;
end
T（k+1,1）=T（k,1）/2+h*Q;
M=M*2;
forj=1:
k
T（k+1,j+1）=T（k+1,j）+（T（k+1,j）-T（k,j））/（4^j-1）;
end
tol=abs（T（k+1,j+1）-T（k,j））;
end
I=T（k+1,k+1）;
step=k;
程序：
[q,
升级会员 