计算方法上机实践一二四文档格式.docx
《计算方法上机实践一二四文档格式.docx》由会员分享，可在线阅读，更多相关《计算方法上机实践一二四文档格式.docx（13页珍藏版）》请在冰点文库上搜索。
%初值
I=zeros（N,1）;
I
（1）=-a*I0+1;
fork=2:
N
I（k）=-a*I（k-1）+1/k;
end
%--------------------------------------------
%[方案II]用递推公式
%I（k-1）=（-I（k）+1/k）/a
II=zeros（N,1）;
ifa>
=N/（N+1）
II（N）=0.5*（1/（（a+1）*（N+1））+1/（a*（N+1）））;
else
II（N）=0.5*（1/（（a+1）*（N+1））+1/N）;
fork=N:
-1:
2
II（k-1）=（-I（k）+1/k）/a;
end
%调用matlab高精度数值积分命令quadl计算以便比较
III=zeros（N,1）;
fork=1:
n=k;
III（k）=quadl（@f,0,1）;
%显示计算结果
%
clc
fprintf（'
\n方案I结果方案II结果精确值'
）
N,
fprintf（'
\nI（%2.0f）%17.7f%17.7f%17.7f'
k,I（k）,II（k）,III（k））
%%把下面的这段保存为f.m
functiony=f（x）%定义函数
globalna%参量n为全局变量
y=x.^n./（a+x）;
%注意:
这里一定要'
点'
运算
return
实验结果：
x1=500000000x2=0
X1=500000000X2=1
a=0.05
方案I结果方案II结果精确值
I
（1）0.84777390.84777390.8477739
I
（2）0.45761130.45761130.4576113
I（3）0.31045280.31045280.3104528
I（4）0.23447740.23447740.2344776
I（5）0.18827610.18827610.1882761
I（6）0.15725290.15725290.1572529
I（7）0.13499450.13499450.1349945
I（8）0.11825030.11825030.1182503
I（9）0.10519860.10519860.1051986
I（10）0.09474010.09474010.0947401
I（11）0.08617210.08617210.0861724
I（12）0.07902470.07902470.0790247
I（13）0.07297180.07297180.0729718
I（14）0.06778000.06778000.0677800
I（15）0.06327770.06327770.0632777
I（16）0.05933610.05933610.0593361
I（17）0.05585670.05585670.0558567
I（18）0.05276270.05276270.0527627
I（19）0.04999340.04999340.0499934
I（20）0.04750030.04767570.0475003
------------------------------------------------------------------------------------------------------
a=15
I
（1）0.03192220.03192220.0319222
I
（2）0.02116730.02116730.0211673
I（3）0.01582450.01582450.0158245
I（4）0.01263260.01263260.0126326
I（5）0.01051120.01051120.0105112
I（6）0.00899930.00899930.0089993
I（7）0.00786740.00786740.0078674
I（8）0.00698830.00698830.0069883
I（9）0.00628620.00628620.0062859
I（10）0.00570640.00570640.0057117
I（11）0.00531360.00531360.0052337
I（12）0.00362890.00362890.0048296
I（13）0.02248960.02248960.0044838
I（14）-0.2659159-0.26591590.0041831
I（15）4.05540504.05540500.0039207
I（16）-60.7685756-60.76857560.0036893
I（17）911.5874579911.58745790.0034837
I（18）-13673.7563129-13673.75631290.0032998
I（19）205106.3973248205106.39732480.0031344
I（20）-3076595.90987240.00307540.0029847
--------------------------------------------------------------------------------------------------------
五实验心得与体会：
第二次实验
实验目的：
1学习非线性方程f（x）=0数值求根命令fzero；
学习用割线法求近似根；
学习绘制隐函数的图像
2学习非线性方程组F（X）=0数值求根命令fsolve。
实验课题：
求解非线性方程（组）的根
采用的具体算法：
functionmysecant
f=inline（'
x-exp（-x）'
）;
a=0.5;
b=0.6;
delta=0.001;
epsilon=0.00000000000001;
max1=50;
[c,err,iter,yc]=secant（f,a,b,delta,epsilon,max1）
%%---------------------------------------------------------
function[c,err,iter,yc]=secant（f,a,b,delta,epsilon,max1）
%%[c,err,iter,yc]=secant（f,a,b,delta,epsilon,max1）
%%输入:
f连续函数
%%a,b迭代初值
%%delta,epsilon容差
%%max1最大迭代次数
%%输出:
c近似根
%%err误差
%%iter迭代次数
%%yc=f（c）
max1
ya=feval（f,a）;
%ya=f（a）
yb=feval（f,b）;
c=a-ya*（b-a）/（yb-ya）;
%割线与x轴交点的横坐标
err=abs（c-b）;
%相邻两次迭代的误差
relerr=err/（abs（c）+eps）;
%相对误差,eps是matlab常数（机器精度）约为1e-16
yc=feval（f,c）;
if（（err<
delta）|（relerr<
delta）|（abs（yc）<
epsilon））%'
|'
是'
或'
break
a=c;
b=b;
iter=k;
functionimplicit_function
globalp%定义全局变量
n=101;
x=linspace（-5,5,101）;
y=zeros（1,n）;
%定义矩阵,初值是零,这是最常用的定义矩阵的方法
y0=-4.6;
%第一次迭代初值
n
p=x（k）;
y（k）=fzero（@fun,y0）;
y0=y（k）;
plot（x,y）%作图
title（'
y^3/（2+0.1*sin（x*y））+x^2-4*x'
）%加个标题
%%------------------------------
functionz=fun（y）%定义函数,这是最常用的定义函数的方式
globalp
x=p;
z=y^3/（2+0.1*sin（x*y））+x^2-4*x;
clf%清图像
x=-1.5:
0.01:
1.5;
%离散化,步长0.01,这也是常用的方法
y1=-（x-3）/2;
%第一个方程求函数值,注意所有运算都是'
运算,和常数相乘等就不用写'
了
y2=（-2*x.^2+5）.^0.5;
%第二个方程求函数值,注意这里自变量与函数是颠倒的
plot（x,y1,x,y2）%作图,注意再把自变量与函数颠倒过来
gridon
axis（[-22-14]）
%限制横坐标与纵坐标的范围,这里要通过不断偿试来得到合适的范围
%通过作图发现在（-1,2）附近有一个根,调用fsolve求更精确的解
clc
X0=[-1,2];
X=fsolve（@myfun,X0）
%------------定义函数组myfun----------------
functionF=myfun（X）
F=[X
（1）-3+2*X
（2）;
2*X
（1）^2+X
（2）^2-5];
iter=5
yc=7.7005e-013
c=0.5671
err=0.0329
iter=6
yc=-4.5519e-015
实验：
X=
-0.82141.9107
心得与体会：
第三次实验
解线性方程组左除命令\的学习
与矩阵计算有关的一些常用函数
病态矩阵试验
OR迭代法收敛速度受松驰因子的影响试验
具体算法：
exp3_1
%【实验二】x=A\b与x=inv（A）*b在耗时方面的区别
%解方程组Ax=b时尽量不要使用x=inv（A）*b,而要使用x=A\b
%下面是二者耗时方面的区别实验（把下面程序拷贝为新的M-文件,有些命令可不关心其含义）
rand（'
state'
0）;
A=gallery（'
randsvd'
200,2e13,2）;
%产生条件数为2e13的200阶的随机矩阵
x=ones（200,1）;
%设精确解为[1,1,...,1]'
b=A*x;
formatlong%用15位小数显示
%%求逆法
tic%启动计时器
x1=inv（A）*b;
time1=toc%关闭计时器并显示耗时
error1=norm（x-x1,inf）%计算最大误差
%%左除法
tic
x2=A\b;
time2=toc
error2=norm（x-x2,inf）
time1=
0.179745959332820
error1=
0.019363403320313
time2=
0.005831327724613
error2=
0.005499554430242
n=13
A的条件数=365015202015855100.000000
Gauss法求解结果pcg法求解结果
1.001.00
1.011.01
0.811.00
2.721.00
-8.321.00
33.461.00
-74.161.00
117.921.00
-119.771.00
80.421.00
-29.091.00
6.000.99
---------------------------------
n=14
A的条件数=4084634689031585800.000000
1.001.01
1.031.00
0.830.99
1.411.00
1.691.00
-6.861.00
26.481.00
-45.461.01
53.341.00
-35.211.00
15.141.00
-1.390.99
实验二取
V=1:
3:
30
c2=1.8063e+015
第四次实验
1学习一维插值命令；
2通过Lagrange插值的Runge现象和Newton插值的Runge现象演示增加对龙格现象的认识；
3学习三次样条的，通过三次样条插值的龙格现象的演示增加对样条插值优点的感性认识。
龙格现象的发生、防止和插值效果的对比
Lagrange插值
样条插值
升级会员 