附：汽轮机课程设计MATLAB计算程序代码：
一、汽轮机流量计算的代码：
%-------------------------------------------------------------------------
%汽轮机流量的计算
P_el=2.25;%汽轮机的的设计功率
deltah_n=200; %汽轮机调节级的比焓降
h0=3224.5;
h_l=h0-deltah_n; %调节级后的比焓
h_e1=2711.4;
h_e2=2637.6;
h_z=2429.74;
h_ed=437;
h_cw=316.6;
h_w2=304.03;
h_w1=171.73;
deltahmac=h0-h_z ; %通流部分的理想比焓降
eta_ri=0.775;
eta_m=0.97;
eta_g=0.97;
eta_h=0.98;
m=1.026
D_0=(3.6*P_el*1000/(deltahmac*eta_g*eta_m)*m)/0.97 %汽轮机的总进气量deltaD_1=0.03*D_0
D_fw=D_0-deltaD_1
deltaD_e1=-D_fw*(-h_ed+h_cw)/(h_e1-h_cw)
deltaD_cw=D_fw*(h_e1-h_ed)/(h_e1-h_cw)
deltaD_e2=deltaD_cw*(h_w2-h_w1)/((h_e2-h_cw)*eta_h)
%调节级
P_i0=D_0*deltah_n/3.6
%第一级组
D_1=D_0-deltaD_1
P_i1=D_1*(h_l-h_e1)/3.6
%第二级组
D_2=D_1-deltaD_e1
P_i2=D_2*(h_e1-h_e2)/3.6
%第三级组
D_3=D_2-deltaD_e2
P_i3=D_3*(h_e2-h_z)/3.6
%整机功率
P_i=P_i0+P_i1+P_i2+P_i3
%机械损失
deltaP_m=P_i*(1-eta_m)
%汽轮机轴端功率
P_a=P_i-deltaP_m
%发电机功率
P_e=P_a*eta_g
%误差计算
erro1=P_el*1000-P_e
erro=(P_el-P_e/1000)/P_el*100
%------------------------------------------------------------------------
%汽耗率
d=D_0*1000/P_e
d_1=D_0*1000/(D_0*(h0-h_z)/3.6-deltaP_m)/eta_g
q=d*(h0-h_ed)
%汽轮机装置的绝对电效率
eta_el=3600/q
%------------------------------------------------------------------------- 二、汽轮机调节级热力计算的代码：
%汽轮机调节级的计算
%-------------------------------------------------------------------------- n=5600;
D_0=11.36;
p_0=2.256;
h_0=3220.17;
v_0=0.1259;
deltah_t=200;
x_a=0.24;
aerfa_1=14;
oum_b=0;
oum_gb=0.03;
oum_b1=0.06;
delta_t=1;
miu_t=0.6;
e_c=0.4
%-------------------------------------------------------------------------- oum_m=oum_b+oum_gb+oum_b1
deltah_n=deltah_t*(1-oum_m)
h_1t=h_0-deltah_n
p_1=input('喷嘴出口压力')
v_1t=input('喷嘴出口比体积')
efthelen=p_1/p_0
c_1t=sqrt(2000*deltah_n)
c_1=0.97*c_1t
k=1.3
sinaerfa1_jia_s1=(2/(k+1))^(1/(k-1))*sqrt((k-1)/(k+1))/(efthelen^(1/k)*sqrt(1-efthelen^((k-1)/ k)))*sin(aerfa_1/180*pi)
s1=asin(sinaerfa1_jia_s1)*180/pi-aerfa_1
u=x_a*c_1
%调节级平均直径
d_m=60*x_a*c_1/(pi*n) %m
%喷管出口面积
A_n=D_0/3.6/(0.648*sqrt(p_0*10^6/v_0))*10000 %cm2
e=0.2% ---------------
l_n=A_n./(e*pi*d_m*100*sin(aerfa_1/180*pi))
deltah_n_kthe=(1-0.97^2)*deltah_n
%--------------------------------------------------------------------------
%第一列动叶热力计算
w_1=sqrt(c_1^2+u^2-2*c_1*u*cos((aerfa_1+s1)/180*pi))
beta_1=asin(c_1*sin((aerfa_1+s1)/180*pi)/w_1)/pi*180
w_2t=w_1
fai_b=-0.07/480*w_2t+0.92+0.7/24
w_2=fai_b.*w_2t
delta_beta_1=input('动叶出口气流角变化，3~5')
beta_2=beta_1-delta_beta_1
c_2=sqrt(w_2.^2+u.^2-2.*w_2.*u.*cos(beta_2/180*pi))
aerfa_2=asin(w_2*sin(beta_2/180*pi)/c_2)/pi*180
h_1=h_1t+deltah_n_kthe %动叶进口比焓
v_1=input('动叶进口比体积，1/5.81')
delta_r=0.5;
delta_t=1.5;
delta=delta_r+delta_t
l_b11=l_n*10+delta_r+delta_t
miu_b=input('查图1-11，0.93')
A_b=D_0/3.6/(miu_b*w_2t/v_1)*10000
l_b1=A_b./(e*pi*d_m*sin(beta_2/180*pi))/10
l_b1-l_b11
pause
deltah_bkthe=w_2t^2/2000*(1-fai_b^2)
h_2=h_1+deltah_bkthe
v_2=input('动叶出口比体积v_2')
%--------------------------------------------------------------------------
%导叶热力计算
deltah_gb=oum_gb*deltah_t
h_1t1=h_2-deltah_gb
p_11=input('导叶出口压力mpa')
v_11=input('导叶出口比体积')
c_1t1=sqrt(2000*deltah_gb+c_2^2)
fai_gb=input('查图1-18,fai_gb')
c_11=fai_gb*c_1t1
delta_aerfa=input('输入deltaaerfa 5~10')
aerfa_11=aerfa_2-delta_aerfa
l_gb1=l_b1+delta
deltaG_gbt=e*0.6*pi*(d_m+l_gb1/1000)*0.001*sqrt(2000*oum_gb*deltah_t)/v_11
G_gb=D_0/3.6-deltaG_gbt
miu_gb=input('查图11-11')
A_gb=G_gb/(miu_gb*c_1t1/v_11)*10000
l_gb=A_gb/(e*pi*d_m*10*sin(aerfa_11/180*pi))
deltah_gb1=c_1t1^2/2000*(1-fai_gb^2)
%导叶出口实际参数
h_1_shiji=h_1t1+deltah_gb1
v_1_shiji=input('由实际出口焓查焓熵图—实际出口比体积')
%--------------------------------------------------------------------------
%第二列动叶热力计算
deltah_b1=oum_b1*deltah_t
h_2t1=h_1_shiji-deltah_b1
p_21=input('动叶出口压力')
v_2t1=input('动叶出口比体积')
w_11=sqrt(c_11^2+u^2-2*c_11*u*cos((aerfa_11)/180*pi))
beta_11=asin(c_11*sin((aerfa_11)/180*pi)/w_11)/pi*180
w_2t1=sqrt(2000*deltah_b1+w_11^2)
fai_b1=input('查图1-18')
w_21=fai_b1*w_2t1
delta_beta2=input('7~18')
beta_21=beta_11-delta_beta2
c_21=sqrt(w_21^2+u^2-2*w_21*u*cos(beta_21/180*pi))
aerfa_21=asin((w_21*sin(beta_21/180*pi))/c_21)*180/pi
deltah_b11=w_2t1^2/2000*(1-fai_b1^2)%动叶损失
deltah_c21=c_21^2/2000
h_21_shiji=h_2t1+deltah_c21
v_21_shiji=input('动叶出口比体积')
l_b21=l_gb+delta
oum_br1=1-(1-oum_b1)*d_m/(d_m-l_b21/1000)
oum_bt1=1-(1-oum_br1)*(1-l_b21/1000)/(1+l_b21/1000)
deltaG_bt1=e*0.6*pi*(d_m+l_b21/1000)*0.001*sqrt(2000*oum_bt1*deltah_t)/v_2t1
miu_b1=input('查图11-11')
A_b1=(D_0/3.6-deltaG_bt1)/(miu_b1*w_2t1/v_2t1)*10000
l_b2=A_b1/(e*pi*d_m*sin(beta_21/180*pi))/10
l_b2-l_b21
%--------------------------------------------------------------------------
%轮周功校核
P_u1_1=u*(c_1*cos((aerfa_1+s1)/180*pi)+c_11*cos(aerfa_11/180*pi)+c_2*cos(aerfa_2/180 *pi)+c_21*cos(aerfa_21/180*pi))/1000
deltah_u=deltah_n_kthe+deltah_bkthe+deltah_gb1+deltah_b11+deltah_c21
P_u1_2=deltah_t-deltah_u
delta_eta=(P_u1_2-P_u1_1)/P_u1_2*100
%--------------------------------------------------------------------------
%轮周效率
eta_u=deltah_u/deltah_t
%--------------------------------------------------------------------------
%级内损失计算
deltaP_f=1.2*(u/100)^3*d_m^2*0.5*(1/v_1+1/v_21_shiji)
deltah_f=3600*deltaP_f/(D_0*1000)
l=(l_n+l_gb1+l_gb+l_b11+l_b1+l_b2+l_b21)/7
deltah_e=2/l*P_u1_2
E_0=deltah_t
x_a1=u/sqrt(2000*deltah_t)
kthe_w=0.55*1/e*(1-e-e_c/2)*x_a1^3
deltah_w=kthe_w*E_0
kthe_s=0.016*1/e*2/d_m*x_a1
deltah_s=kthe_s*E_0
deltah_e1=deltah_w+deltah_s
deltahh=deltah_t-deltah_n_kthe-deltah_bkthe-deltah_gb1-deltah_b11-deltah_e-deltah_c21 deltah_t1=(deltaG_gbt+deltaG_bt1)/(D_0/3.6)*deltahh
%--------------------------------------------------------------------------
%级的内功率
deltah_i=deltahh-deltah_f-deltah_t1-deltah_e1
P_i=D_0/3.6*deltah_i
%--------------------------------------------------------------------------
%级的内效率
eta_i=deltah_i/E_0
%--------------------------------------------------------------------------
三、汽轮机压力级热力计算的代码：
d_m=input('级的平均直径');
c_0=input('流入该级的蒸汽初速');
D=input('流量');
p_0=input('级前蒸汽压力');
x_0=input('干度');
deltah_t=input('级的理想比焓');
oum_m=input('级的平均反动度');
oum_t=input('顶部反动度');
aerfa_1=input('喷管出气角');
miu_1=input('余速利用系数');
%喷管的理想滞止比焓降
deltah_c0=c_0^2/2000
deltah_n_xing=(1-oum_m)*(deltah_t+deltah_c0)
deltah_t_xing=deltah_t+deltah_c0
deltah_b=oum_m*deltah_t_xing
%喷管热力计算
p_0_xing=input('喷管前的滞止压力');
h_0_xing=input('滞止比焓');
v_0_xing=input('滞止比体积');
h_0=input('喷管前比焓');
p_1=input('喷管后压力');
v_1t=input('理想比体积');
h_1t=h_0-deltah_t
k=1.035+0.1*x_0
efthlen_cr=(2/(k+1))^(k/(k-1))
efthlen_n=p_1/p_0_xing
pause
p_cr=p_0_xing*efthlen_cr
h_cr=input('临界焓');
v_cr=input('临界比体积');
c_cr=sqrt(2000*(h_0_xing-h_cr))
c_1t=sqrt(2000*deltah_n_xing)
fai=0.97;
c_1=fai*c_1t
aerfa_jia_delta=asin(sin(aerfa_1/180*pi)*c_cr*v_1t/(c_1t*v_cr))/pi*180
miu_p=input('0.8~0.7');
d_p=input('汽封齿的平均直径');
delta_p=input('汽封间隙');
z_p=input('齿数');
deltaD_p=miu_p*pi*d_p*delta_p*sqrt(2*deltah_n_xing*1000)/(v_1t*sqrt(z_p+1/z_p)) D_n=D/3.6-deltaD_p
A_n=D_n/(0.648*sqrt(p_0_xing*10^6/v_0_xing))
l_n=A_n/(1*pi*d_m*sin(aerfa_1/180*pi))*1000
deltah_nkthe=(1-fai^2)*deltah_n_xing
h_1=h_1t+deltah_nkthe
v_1=input('喷管出口实际比体积');
%---------------------------------------------
%动叶栅热力计算
n=3000
u=pi*d_m*n/60
w_1=sqrt(c_1^2+u^2-2*c_1*u*cos(aerfa_jia_delta/180*pi))
beta_1=atan(c_1*sin(aerfa_jia_delta/180*pi)/(c_1*cos(aerfa_jia_delta/180*pi)-u))/pi*1 80
deltah_w1=w_1^2/2000
w_2t=sqrt(2000*(deltah_b+deltah_w1))
faifai=input('速度系数');
w_2=faifai*w_2t
beta_2=input('beta2');
c_2=sqrt(w_2^2+u^2-2*w_2*u*cos(beta_2/180*pi))
aerfa_2=asin(w_2*sin(beta_2*pi/180)/c_2)/pi*180
deltah_bkthe=w_2^2/2000*(1-faifai^2)
h_2t=h_1-deltah_b
v_2t=input('动叶出口理想比体积');
p_2=input('动叶出口压力');
efthlen_b=p_2/p_1
pause
h_2=h_2t+deltah_bkthe
v_2=input('动叶出口实际比体积');
x_2=input('动叶出口实际干度');
l_b1=l_n+1+2
deltaD_t=1*0.6*pi*(1.442+0.075)*0.002/v_2t*sqrt(oum_t*deltah_t_xing*1000)
miu_b=0.935;
A_b=(D_n-deltaD_t)/(miu_b*w_2t/v_2t)
d_b=input('输入直径db')
l_b=A_b/(pi*d_b*sin(beta_2*pi/180))*1000
deltah_c2=c_2^2/2000
%-------------------------------------------------------------------
%轮周效率
P_u=deltah_t_xing-deltah_nkthe-deltah_bkthe-deltah_c2
P_u1=u*(c_1*cos(aerfa_jia_delta/180*pi)+c_2*cos(aerfa_2*pi/180))/1000
deltah_u=P_u;
eta_u=deltah_u/(deltah_t_xing-miu_1*deltah_c2)
%级内各项损失计算
a=input('输入a的值');
deltah_l=a/l_n*deltah_u
theta=input('输入theta');
k_1=input('输入k1');
d=(d_m+d_b)/2;
deltap_f=k_1*(u/100)^3*d^2*1/theta
deltah_f=deltap_f/D*3.6
deltah_e=input('deltahe');
deltah_u1=deltah_t_xing-deltah_nkthe-deltah_c2-deltah_e
deltah_p=deltaD_p/(D/3.6)*deltah_u1
deltah_t=deltaD_t/(D/3.6)*deltah_u1
deltah_f=input('deltah_f');
deltah_i1=deltah_t_xing-deltah_nkthe-deltah_c2-deltah_l-deltah_p-deltah_t-deltah_f deltah_x=(1-(x_0+x_2)/2)*deltah_i1
deltah_i=deltah_i1-deltah_x
P_i=D/3.6*deltah_i
eta_ri=deltah_i/(deltah_t_xing-miu_1*deltah_c2)
%----------------------------------------------------------------
6-1 试验目的计算特征值，实现算法 试验容：随机产生一个10阶整数矩阵，各数均在-5和5之间。 （1） 用MATLAB 函数“eig ”求矩阵全部特征值。 （2） 用幂法求A 的主特征值及对应的特征向量。 （3） 用基本QR 算法求全部特征值（可用MATLAB 函数“qr ”实现矩阵的QR 分解）。 原理 幂法：设矩阵A 的特征值为12n ||>||||λλλ≥???≥并设A 有完全的特征向量系12,,,n χχχ???(它们线性无关)，则对任意一个非零向量0n V R ∈所构造的向量序列1k k V AV -=有11()lim ()k j k k j V V λ→∞ -=， 其中()k j V 表示向量的第j 个分量。 为避免逐次迭代向量k V 不为零的分量变得很大（1||1λ>时）或很小（1||1λ<时），将每一步的k V 按其模最大的元素进行归一化。具体过程如下： 选择初始向量0V ，令1max(),,,1k k k k k k k V m V U V AU k m +===≥，当k 充分大时1111,max()max() k k U V χλχ+≈ ≈。 QR 法求全部特征值： 111 11222 111 ,1,2,3,k k k k k A A Q R R Q A Q R k R Q A Q R +++==????==??=???? ??????==?? 由于此题的矩阵是10阶的，上述算法计算时间过长，考虑采用改进算法——移位加速。迭 代格式如下： 1 k k k k k k k k A q I Q R A R Q q I +-=?? =+? 计算k A 右下角的二阶矩阵() () 1,1 1,() (),1 ,k k n n n n k k n n n n a a a a ----?? ? ??? 的特征值()()1,k k n n λλ-，当()()1,k k n n λλ-为实数时，选k q 为()()1,k k n n λλ-中最接近(),k n n a 的。 程序
1 单选(2分) 利用MATLAB求解科学计算问题的优势是（）。 得分/总分 ? A. 算法最优 ? B. 不需要编写程序 ? C. 程序执行效率高 ? D. 编程效率高 正确答案：D你没选择任何选项 2 单选(2分) 在MATLAB命令行窗口输入命令时，可使用续行符，其写法是（）。 得分/总分 ? A. 省略号（…） ? B. 分号（;） ? C. 三个小数点（…） ? D. 百分号（%） 正确答案：C你没选择任何选项 3
下列语句执行后，D的值为（）。 1.A=[1:3;4:6]; 2.D=sub2ind(size(A),[1,1],[2,3]) 得分/总分 ? A. 3 6 ? B. 2 5 ? C. 3 5 ? D. 4 5 正确答案：C你没选择任何选项 4 单选(2分) ceil(-2.1)+floor(-2.1)+fix(-2.1)的结果为（）。 得分/总分 ? A. -7 ? B. -6 ? C. -5 ? D. -9 正确答案：A你没选择任何选项 5
下列语句执行后，x的值是（）。 1.log=1:5; 2.x=log(1) 得分/总分 ? A. ? B. 1 ? C. 数学常数e ? D. 报错 正确答案：B你没选择任何选项 6 单选(2分) 下列语句执行后，c的值是（）。 1.ch=['abcdef';'123456']; 2.c=char(ch(2,4)-1) 得分/总分 ? A. '4' ? B. 4 ? C. '3' ? D. 3
7 单选(2分) 产生和A同样大小的全0矩阵的函数是（）。 得分/总分 ? A. zero(size(A)) ? B. zeros(size(A)) ? C. size(zero(A)) ? D. size(zeros(A)) 正确答案：B你没选择任何选项 8 单选(2分) 语句x=speye(5)==eye(5)执行后，则下列说法中正确的是（）。 得分/总分 ? A. x是5阶全1矩阵，且采用稀疏存储方式 ? B. x是5阶全1矩阵，且采用完全存储方式 ? C. x是5阶单位矩阵，且采用稀疏存储方式 ? D. x是5阶单位矩阵，且采用完全存储方式
%求四阶线性方程组的MA TLAB程序 clear Ab=[0.001 2 1 5 1; 3 - 4 0.1 -2 2; 2 -1 2 0.01 3; 1.1 6 2.3 9 4];%增广矩阵 num=[1 2 3 4];%未知量x的对应序号 for i=1:3 A=abs(Ab(i:4,i:4));%系数矩阵取绝对值 [r,c]=find(A==max(A(:))); r=r+i-1;%最大值对应行号 c=c+i-1;%最大值对应列号 q=Ab(r,:),Ab(r,:)=Ab(i,:),Ab(i,:)=q;%行变换 w=Ab(:,c),Ab(:,c)=Ab(:,i),Ab(:,i)=w;%列变换 n=num(i),num(i)=num(c),num(c)=n;%列变换引起未知量x次序变化for j=i:3 Ab(j+1,:)=-Ab(j+1,i)*Ab(i,:)/Ab(i,i)+Ab(j+1,:);%消去过程 end end %最后得到系数矩阵为上三角矩阵 %回代算法求解上三角形方程组 x(4)=Ab(4,5)/Ab(4,4); x(3)=(Ab(3,5)-Ab(3,4)*x(4))/Ab(3,3); x(2)=(Ab(2,5)-Ab(2,3)*x(3)-Ab(2,4)*x(4))/Ab(2,2); x(1)=(Ab(1,5)-Ab(1,2)*x(2)-Ab(1,3)*x(3)-Ab(1,4)*x(4))/Ab(1,1); for s=1:4 fprintf('未知量x%g =%g\n',num(s),x(s)) end %验证如下 %A=[0.001 2 1 5 1; 3 -4 0.1 -2 2;2 -1 2 0.01 3; 1.1 6 2.3 9 4]; %b=[1 2 3 4]'; %x=A\b; %x1= 1.0308 %x2= 0.3144 %x3= 0.6267 %x4= -0.0513
一、Maple V 系统 Maple V是由Waterloo大学开发的数学系统软件，它不但具有精确的数值处理功能，而且具有无以伦比的符号计算功能。Maple V的符号计算能力还是MathCAD和MATLAB等软件的符号处理的核心。Maple提供了2000余种数学函数，涉及范围包括：普通数学、高等数学、线性代数、数论、离散数学、图形学。它还提供了一套内置的编程语言，用户可以开发自己的应用程序，而且Maple自身的2000多种函数，基本上是用此语言开发的。 Maple采用字符行输入方式，输入时需要按照规定的格式输入，虽然与一般常见的数学格式不同，但灵活方便，也很容易理解。输出则可以选择字符方式和图形方式，产生的图形结果可以很方便地剪贴到Windows应用程序内。 二、MATLAB 系统 MATLAB原是矩阵实验室（Matrix Laboratory）在70年代用来提供Linpack和Eispac k软件包的接口程序，采用C语言编写。从80年代出现3.0的DOS版本，逐渐成为科技计算、视图交互系统和程序语言。MATLAB可以运行在十几个操作平台上，比较常见的有基于W indows 9X/NT、OS/2、Macintosh、Sun、Unix、Linux等平台的系统。 MATLAB程序主要由主程序和各种工具包组成，其中主程序包含数百个内部核心函数，工具包则包括复杂系统仿真、信号处理工具包、系统识别工具包、优化工具包、神经网络工具包、控制系统工具包、μ分析和综合工具包、样条工具包、符号数学工具包、图像处理工具包、统计工具包等。而且5.x版本还包含一套几十个的PDF文件，从MATLAB的使用入门到其他专题应用均有详细的介绍。 MATLAB是数值计算的先锋，它以矩阵作为基本数据单位，在应用线性代数、数理统计、自动控制、数字信号处理、动态系统仿真方面已经成为首选工具，同时也是科研工作人员和大学生、研究生进行科学研究的得力工具。MATLAB在输入方面也很方便，可以使用内部的E ditor或者其他任何字符处理器，同时它还可以与Word6.0/7.0结合在一起，在Word的页面里直接调用MATLAB的大部分功能，使Word具有特殊的计算能力。 三、MathCAD 系统 MathCAD是美国Mathsoft公司推出的一个交互式的数学系统软件。从早期的DOS下的1. 0和Windows下的4.0版本，到今日的8.0版本，功能也从简单的数值计算，直至引用Map le强大的符号计算能力，使得它发生了一个质的飞跃。 MathCAD是集文本编辑、数学计算、程序编辑和仿真于一体的软件。MathCAD7.0 Profe ssional（专业版）运行在Win9X/NT下，它的主要特点是输入格式与人们习惯的数学书写格式很近似，采用WYSWYG（所见所得）界面，特别适合一般无须进行复杂编程或要求比较特殊的计算。MathCAD 7.0 Professional 还带有一个程序编辑器，对于一般比较短小，或者要求计算速度比较低时，采用它也是可以的。这个程序编辑器的优点是语法特别简单。 MathCAD可以看作是一个功能强大的计算器，没有很复杂的规则；同时它也可以和Wor d、Lotus、WPS2000等字处理软件很好地配合使用，可以把它当作一个出色的全屏幕数学公式编辑器。 四、Mathematica 系统 Mathematica是由美国物理学家Stephen Wolfram领导的Wolfram Research开发的数学系统软件。它拥有强大的数值计算和符号计算能力，在这一方面与Maple类似，但它的符
第四讲绘图功能
作为一个功能强大的工具软件，Matlab 具有很强的图形处理功能，提供了大量的二维、三维图形函数。由于系统采用面向对象的技术和丰富的矩阵运算，所以在图形处理方面即常方便又高效。
4.1 二维图形 一、plot函数 函数格式：plot(x,y)其中x和y为坐标向量函数功能：以向量x、y为轴，绘制曲线。【例1】在区间0≤X≤2 内，绘制正弦曲线Y=SIN（X），其程序为： x=0:pi/100:2*pi; y=sin(x); plot(x,y)
一、plot函数 【例2】同时绘制正、余弦两条曲线Y1=SIN（X）和Y2=COS（X），其程序为： x=0:pi/100:2*pi; y1=sin(x); y2=cos(x); plot(x,y1,x,y2) plot函数还可以为plot(x,y1,x,y2，x,y3，…)形式，其功能是以公共向量x为X轴，分别以y1，y2，y3，…为Y轴，在同一幅图内绘制出多条曲线。
一、plot函数 （一）线型与颜色 格式：plot(x,y1,’cs’,...) 其中c表示颜色，s表示线型。 【例3】用不同线型和颜色重新绘制例4.2图形，其程序为：x=0:pi/100:2*pi; y1=sin(x); y2=cos(x); plot(x,y1,'go',x,y2,'b-.') 其中参数'go'和'b-.'表示图形的颜色和线型。g表示绿色，o表示图形线型为圆圈；b表示蓝色，-.表示图形线型为点划线。
一、plot函数 （二）图形标记 在绘制图形的同时，可以对图形加上一些说明，如图形名称、图形某一部分的含义、坐标说明等，将这些操作称为添加图形标记。 title(‘加图形标题'); xlabel('加X轴标记'); ylabel('加Y轴标记'); text(X,Y,'添加文本');
练习题 1.求函数在指定点的数值导数 x=sym('x'); >> y=[x x.^2 x.^3;1 2*x 3*x.^2;0 2 6*x]; >> x=1; >> eval(diff(y)) ans = 1 2 3 0 2 6 0 0 6 >> x=2; >> eval(diff(y)) ans = 1 4 12 0 2 12 0 0 6 >> x=3; >> eval(diff(y)) ans = 1 6 27 0 2 18 0 0 6 2.求下列函数导数 （1） x=sym('x'); >> y=x^10+10^x+(log(10))/log(x); >> diff(y) ans = 10*x^9+10^x*log(10)-2592480341699211/1125899906842624/log(x)^2/x （2） x=sym('x');
>> y=log(1+x); >> x=1; >> eval(diff(y,2)) %在x=1的条件下对y表达式求两次导数后导函数的值 ans = -0.2500 3.用数值方法求下列积分 首先先讲一下trapz的用法，如下题 t=0:0.001:1; >> y=t; >> trapz(t,y) ans = 0.5000 （1） >> x=1:0.01:5; >> y=(x.^2).*sqrt(2*x.^2+3); >> trapz(x,y) ans = 232.8066 (2) x=pi/4:0.01:pi/3; >> y=x./(sin(x).^2); >> trapz(x,y) ans = 0.3810 第三题拟合曲线题 x=[0:0.1:1]; >> y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2]; >> a=polyfit(x,y,2); >> x=[0.05:0.2:1.05]; >> y=a(3)+a(2)*x+a(1)*x.^2 %注意x要在y前先赋值，不然y不会运行为最新的x对呀的y值 y =
《Matlab与科学计算》作业 第一章MATLAB环境 1、MATLAB通用操作界面窗口包括哪些？命令窗口、历史命令窗口、当前目录窗口、工作空间窗口各有哪些功能？ 答：MATLAB通用操作界面窗口包括：命令窗口、历史命令窗口、当前目录浏览器窗口、工作空间窗口、变量编辑器窗口、M文件编辑/调试器窗口、程序性能剖析窗口、MATLAB帮助。 命令窗口是MATLAB命令操作的最主要窗口，可以把命令窗口当做高级的“草稿纸”。在命令窗口中可以输入各种MATLAB的命令、函数和表达式，并显示除图形外的所有运算结果。 历史命令窗口用来记录并显示已经运行过的命令、函数和表达式，并允许用户对它们进行选择、复制和重运行，用户可以方便地输入和修改命令，选择多行命令以产生M文件。 当前目录窗口用来设置当前目录，可以随时显示当前目录下的M、MKL等文件的信息，扬文件类型、文件名、最后个修改时间和文件的说明信息等，并可以复制、编辑和运行M文件及装载MAT数据文件。 工作空间窗口用来显示所有MATLAB工作空间中的变量名、数据结构、类型、大小和字节数。 2、熟悉课本中表格1.4、1.5、1.6、1.7、1.8的内容。 3、如何生成数据文件？如何把数据文件中的相关内容输入到工作空间中，用实例进行操作。 生成数据文件：
把数据文件中的相关内容输入到工作空间中： 结果： 4、在工作空间中可以通过哪些命令管理变量，写出每种语法的具体操作过程。答：（1）把工作空间中的数据存放到MAT数据文件。 语法：save filename 变量1 变量2 ……参数。 （2）从数据文件中取出变量存放到工作空间。 语法：load filename 变量1 变量2 ……。
盐城师范学院《MATLAB与科学计算》期末论文 2016－2017学年度第一学期 用MATLAB解决解析几何的图形问题 学生姓名吴梦成 学院数学与统计学院 专业信息与计算科学 班级数15(5)信计 学号 15213542
用MATLAB 解决解析几何的图形问题 摘 要 将 MATLAB 的图形和动画功能都用于解析几何教学，可使教学形象生动。以图形问题为例，详细给出了实例的程序编写和动画实现过程 。在解析几何教学中有一定的应用价值。 【关键词】： MATLAB ； 解析几何 ；图形 ； 动 画；编程 1 引 言 在解析几何的教学中，使用传统的教学方法。许多曲线及曲面的形成过程与变换过程只通过传统的教师讲授静态图示就很难形象生动地表示出来 。在解析几何教学中使用MATLAB 软件辅助教学，不仅可以很容易绘制出复杂的立体图形。把曲线、曲面的形成和变化过程准确地模拟出来 ，而且还能够对它们进行翻转 、旋转 ，甚 至还能够轻而易举地实现图形的动画效果!这对提高教学效率和培养学生的空间想象能力可起到事半功倍的效果。下面结合实例从几个方面说明MATLAB 在解析几何画图方面的应用。 2 利用 MATLAB 绘制三维曲线 在空间解析几何中，各种曲线和曲面方程的建立都离不开图形 ，而空间曲线和曲面图形既难画又费时。借助MATLAB 的绘图功能 ，可以快捷 、 准确地绘出图形，使教学变得形象 、生动 。有利于学生观察三维空间图形的形状 ， 掌握图形的性质 。 一 般地 ，MATLAB 可用plot3，ezplot3，comet3等函数来各种三维曲线 。 例如画螺旋曲线的图形，其参数方程设为 ：t at cos x =，t b sin t y -=，ct =z 。使用 plot3语句画螺旋曲线图形的方法如下( 设a =2 ,b=4,c=3)： );*3),sin(*.*4),cos(*.*2(3;*10:50/:0t t t t t plot pi pi t -= MATLAB 用两条简单的语句就可以画出螺旋 曲线(图1)，但上述方法是静态的 ，为了体
第一章插值方法 1.1Lagrange插值 1.2逐步插值 1.3分段三次Hermite插值 1.4分段三次样条插值 第二章数值积分 2.1 Simpson公式 2.2 变步长梯形法 2.3 Romberg加速算法 2.4 三点Gauss公式 第三章常微分方程德差分方法 3.1 改进的Euler方法 3.2 四阶Runge-Kutta方法 3.3 二阶Adams预报校正系统 3.4 改进的四阶Adams预报校正系统 第四章方程求根 4.1 二分法 4.2 开方法 4.3 Newton下山法 4.4 快速弦截法 第五章线性方程组的迭代法 5.1 Jacobi迭代 5.2 Gauss-Seidel迭代 5.3 超松弛迭代 5.4 对称超松弛迭代 第六章线性方程组的直接法 6.1 追赶法 6.2 Cholesky方法 6.3 矩阵分解方法 6.4 Gauss列主元消去法
第一章插值方法 1.1Lagrange插值 计算Lagrange插值多项式在x=x0处的值. MATLAB文件：（文件名：Lagrange_eval.m）function [y0,N]= Lagrange_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Lagrange插值多项式在x0处的值 %N是Lagrange插值函数的权系数 m=length(X); N=zeros(m,1); y0=0; for i=1:m N(i)=1; for j=1:m if j~=i; N(i)=N(i)*(x0-X(j))/(X(i)-X(j)); end end y0=y0+Y(i)*N(i); end 用法》X=[…];Y=[…]; 》x0= ; 》[y0,N]= Lagrange_eval(X,Y,x0) 1.2逐步插值 计算逐步插值多项式在x=x0处的值. MATLAB文件：（文件名：Neville_eval.m）function y0=Neville_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Neville逐步插值多项式在x0处的值 m=length(X); P=zeros(m,1); P1=zeros(m,1); P=Y; for i=1:m P1=P; k=1; for j=i+1:m k=k+1;
单元测验已完成成绩：分 1 【单选题】 MATLAB一词来自（）的缩写。 A、 Mathematica Laboratory B、 Matrix Laboratory C、 MathWorks Lab D、 Matrices Lab 我的答案：B得分：分 2 【单选题】 下列选项中能反应MATLAB特点的是（）。 A、 算法最优 B、 不需要写程序 C、 程序执行效率高 D、 编程效率高
我的答案：D得分：分 单元测验已完成成绩：分 1 【单选题】 当在命令行窗口执行命令时，如果不想立即在命令行窗口中输出结果，可以在命令后加上（）。 A、 冒号（:） B、 逗号（,） C、 分号（;） D、 百分号（%） 我的答案：C得分：分 2 【单选题】 fix(264/100)+mod(264,10)*10的值是（）。 A、 86 B、 62 C、 423 D、
42 我的答案：D得分：分 3 【单选题】 在命令行窗口输入下列命令后，x的值是（）。 >> clear >> x=i*j A、 不确定 B、 -1 C、 1 D、 i*j 我的答案：B得分：分 4 【单选题】 使用语句x=linspace(0,pi,6)生成的是（）个元素的向量。 A、 8 B、 7 C、 6
D、 5 我的答案：C得分：分 5 【单选题】 ceil的结果为（）。 A、 -2 B、 -3 C、 1 D、 2 我的答案：A得分：分 6 【单选题】 eval('sqrt(4)+2')的值是（）。 A、 sqrt(4)+2 B、 4 C、 2 D、
2+2 我的答案：B得分：分 7 【单选题】 已知a为3×5矩阵，则执行完a(:,[2,4])=[]后（）。 A、 a变成行向量 B、 a变为3行2列 C、 a变为3行3列 D、 a变为2行3列 我的答案：C得分：分 8 【单选题】 在命令行窗口输入以下命令 >> A=[1:3;4:6]; >> D=sub2ind(size(A),[1,1],[2,3]) D的值为（）。 A、 3 6 B、 2 5 C、 4 5
常用计算方法 1．超越方程的求解 一超越方程为 x (2ln x – 3) -100 = 0 求超越方程的解。 [算法]方法一：用迭代算法。将方程改为 01002ln()3 x x =- 其中x 0是一个初始值，由此计算终值x 。取最大误差为e = 10-4，当| x - x 0| > e 时，就用x 的值换成x 0的值，重新进行计算；否则| x - x 0| < e 为止。 [程序]P1_1abs.m 如下。 %超越方程的迭代算法 clear %清除变量 x0=30; %初始值 xx=[]; %空向量 while 1 %无限循环 x=100/(2*log(x0)-3); %迭代运算 xx=[xx,x]; %连接结果 if length(xx)>1000,break ,end %如果项数太多则退出循环（暗示发散） if abs(x0-x)<1e-4,break ,end %当精度足够高时退出循环 x0=x; %替换初值 end %结束循环 figure %创建图形窗口 plot(xx,'.-','LineWidth',2,'MarkerSize',12)%画迭代线'.-'表示每个点用.来表示，再用线连接 grid on %加网格 fs=16; %字体大小 title('超越方程的迭代折线','fontsize',fs)%标题 xlabel('\itn','fontsize',fs) %x 标签 ylabel('\itx','fontsize',fs) %y 标签 text(length(xx),xx(end),num2str(xx(end)),'fontsize',fs)%显示结果 [图示]用下标作为自变量画迭代的折线。如P0_20_1图所示，当最大误差为10-4时，需要迭代19次才能达到精度，超越方程的解为27.539。 [算法]方法二：用求零函数和求解函数。将方程改为函数 100()2ln()3f x x x =-- MATLAB 求零函数为fzero ，fzero 函数的格式之一是 x = fzero(f,x0) 其中，f 表示求解的函数文件，x0是估计值。fzero 函数的格式之二是 x = fzero(f,[x1,x2])
Matlab 与科学计算考试样题（客观题） 1 下面的MATLAB 语句中正确的有： a) 2a ＝pi 。 b) record_1=3+4i c) a=2.0, d) c=1+6j 2. 已知水的黏度随温度的变化公式如下，其中a=0.03368，b=0.000221，计算温度t 为20，30，40度时的粘度分别是： 2 1at bt μμ=++0μ为0℃水的黏度，值为31.78510-?；a 、b 为常数，分别为0.03368、0.000221。 3. 请补充语句以画出如图所示的图形： [x,y]=meshgrid(-2:0.1:2, -2:0.1:2)。 Z=x.*exp(-x.^2-y.^2)。 。 a) Plot3(x,y,Z) b) plot3(x,y,Z) c) mesh(x,y,Z) d) plot3(x,y,z) 2 a) 0.4900 1.2501 0.8560 b) 0.8560 1.2501 0.4900 c) -0.6341 3.8189 -3.7749 d) 3.8189 -3.7749 2.8533 解释说明：
>> x=0.5:0.5:3.0。 >> y=[1.75,2.45,3.81,4.80,8.00,8.60]。 >> a=polyfit(x,y,2) a = 0.4900 1.2501 0.8560 >> x1=[0.5:0.25:3.0]。 >> y1=a(1)*x1.^2+a(2)*x1+a(3) >> plot(x,y,'*') >> hold on >> plot(x1,y1,'--r') 5. 求方程在 x=0.5附近的根. 21 x x += a) 0.6180 b) -1.1719e-25 c) －1 d) -1.6180 6. 用Newton-Cotes方法计算如下积分 1 5 x? （a）133.6625 (b)23.8600 (c) 87.9027 (d) -1.6180 7. y=ln(1+x)，求x=1时y" a) -0.25 b) 0.5 c) -0.6137 d) -1.6137 8.某公司用3台轧机来生产规格相同的铝合金薄板。取样测量薄板的 厚度，精确至‰厘M。得结果如下： 轧机1：0.236 0.238 0.248 0.245 0.243 轧机2：0.257 0.253 0.255 0.254 0.261 轧机3：0.258 0.264 0.259 0.267 0.262 计算方差分析结果，并判定各台轧机所生产的薄板的厚度有无显著的差异？ a) p＝1.3431e-005，没有显著差异。
1、Newdon迭代法求解非线性方程 function [x k t]=NewdonToEquation(f,df,x0,eps) %牛顿迭代法解线性方程 %[x k t]=NewdonToEquation(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:原函数，定义为内联函数 ?:函数的倒数，定义为内联函数 %x0:初始值 %eps:误差限 % %应用举例： %f=inline('x^3+4*x^2-10'); ?=inline('3*x^2+8*x'); %x=NewdonToEquation(f,df,1,0.5e-6) %[x k]=NewdonToEquation(f,df,1,0.5e-6) %[x k t]=NewdonToEquation(f,df,1,0.5e-6) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquation(f,df,1) if nargin==3 eps="0".5e-6; end tic; k=0; while 1 x="x0-f"(x0)./df(x0); k="k"+1; if abs(x-x0) < eps || k >30 break; end x0=x; end t=toc; if k >= 30 disp('迭代次数太多。'); x="0"; t="0"; end
2、Newdon迭代法求解非线性方程组 function y="NewdonF"(x) %牛顿迭代法解非线性方程组的测试函数 %定义是必须定义为列向量 y(1,1)=x(1).^2-10*x(1)+x(2).^2+8; y(2,1)=x(1).*x(2).^2+x(1)-10*x(2)+8; return; function y="NewdonDF"(x) %牛顿迭代法解非线性方程组的测试函数的导数 y(1,1)=2*x(1)-10; y(1,2)=2*x(2); y(2,1)=x(2).^+1; y(2,2)=2*x(1).*x(2)-10; return; 以上两个函数仅供下面程序的测试 function [x k t]=NewdonToEquations(f,df,x0,eps) %牛顿迭代法解非线性方程组 %[x k t]=NewdonToEquations(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:方程组（事先定义） ?:方程组的导数（事先定义） %x0:初始值 %eps:误差限 % %说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义% 另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示 % %应用举例： %x0=[0,0];eps=0.5e-6; %x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)
第8章 常微分方程初值问题数值解法 8.1 引言 8.2 欧拉方法 8.3 龙格-库塔方法 8.4 单步法的收敛性与稳定性 8.5 线性多步法
8.1 引 言 考虑一阶常微分方程的初值问题 00(,),[,],(). y f x y x a b y x y '=∈=（1.1） （1.2） 如果存在实数 ，使得 121212(,)(,).,R f x y f x y L y y y y -≤-?∈（1.3） 则称 关于 满足李普希茨(Lipschitz ）条件， 称为 的李普希茨常数（简称Lips.常数）. 0>L f y L f (参阅教材386页)
计算方法及MATLAB 实现 所谓数值解法，就是寻求解 在一系列离散节点 )(x y <<<<<+121n n x x x x 上的近似值 . ,,,,,121+n n y y y y 相邻两个节点的间距 称为步长. n n n x x h -=+1 如不特别说明，总是假定 为定数， ),2,1( ==i h h i 这时节点为 . ) ,2,1,0(0 =+=i nh x x n 初值问题(1.1),(1.2)的数值解法的基本特点是采取 “步进式”. 即求解过程顺着节点排列的次序一步一步地向前推进. 00(,),[,], (). y f x y x a b y x y '=∈=
描述这类算法，只要给出用已知信息 ,,,21--n n n y y y 计算 的递推公式. 1+n y 一类是计算 时只用到前一点的值 ，称为单步法. 1+n y n y 另一类是用到 前面 点的值 ， 1+n y k 11,,,+--k n n n y y y 称为 步法. k 其次，要研究公式的局部截断误差和阶，数值解 与 精确解 的误差估计及收敛性，还有递推公式的计算 稳定性等问题. n y )(n x y 首先对方程 离散化，建立求数值解的递推 公式. ),(y x f y ='
单元测验已完成成绩：100.0分 1 【单选题】 MATLAB一词来自（）的缩写。 ?A、 Mathematica Laboratory ?B、 Matrix Laboratory ?C、 MathWorks Lab ?D、 Matrices Lab 我的答案：B得分：50.0分 2 【单选题】 下列选项中能反应MATLAB特点的是（）。 ?A、 算法最优 ?B、 不需要写程序 ?C、 程序执行效率高 . .
?D、 编程效率高 我的答案：D得分：50.0分 单元测验已完成成绩：96.4分 1 【单选题】 当在命令行窗口执行命令时，如果不想立即在命令行窗口中输出结果，可以在命令后加上（）。 ?A、 冒号（:） ?B、 逗号（,） ?C、 分号（;） ?D、 百分号（%） 我的答案：C得分：7.1分 2 【单选题】 fix(264/100)+mod(264,10)*10的值是（）。 ?A、 86 . .
?B、 62 ?C、 423 ?D、 42 我的答案：D得分：7.1分 3 【单选题】 在命令行窗口输入下列命令后，x的值是（）。 >> clear >> x=i*j ?A、 不确定 ?B、 -1 ?C、 1 ?D、 i*j 我的答案：B得分：7.1分 . .
4 【单选题】 使用语句x=linspace(0,pi,6)生成的是（）个元素的向量。 ?A、 8 ?B、 7 ?C、 6 ?D、 5 我的答案：C得分：7.1分 5 【单选题】 ceil(-2.1)的结果为（）。 ?A、 -2 ?B、 -3 ?C、 1 . .
?D、 2 我的答案：A得分：7.1分 6 【单选题】 eval('sqrt(4)+2')的值是（）。 ?A、 sqrt(4)+2 ?B、 4 ?C、 2 ?D、 2+2 我的答案：B得分：7.1分 7 【单选题】 已知a为3×5矩阵，则执行完a(:,[2,4])=[]后（）。 ?A、 a变成行向量 ?B、 a变为3行2列 . .
用MATLAB实现结构可靠度计算 口徐华…朝泽刚‘u刘勇‘21 。 (【l】中国地质大学(武汉工程学院湖北?武汉430074; 12】河海大学土木工程学院江苏?南京210098 摘要:Matlab提供了各种矩阵的运算和操作,其中包含结构可靠度计算中常用的各种数值计算方法工具箱,本文从基本原理和相关算例分析两方面,阐述利用Matlab,编制了计算结构可靠度Matlab程.序,使得Matlab-语言在可靠度计算中得到应用。 关键词:结构可靠度Matlab软件最优化法 中图分类号:TP39文献标识码:A文章编号:1007-3973(200902-095-Ol 1结构可靠度的计算方法 当川概率描述结构的可靠性时,计算结构可靠度就是计算结构在规定时问内、规定条件F结构能够完成预定功能的概率。 从简单到复杂或精确稃度的不同,先后提出的可靠度计算方法有一次二阶矩方法、二次二阶矩方法、蒙特卡洛方法以及其他方法。一次■阶矩方法又分为。I-心点法和验算点法,其中验算点法足H前可靠度分析最常川的方法。 2最优化方法计算可靠度指标数学模型 由结构111n个任意分布的独立随机变量一,x:…以表示的结构极限状态方程为:Z=g(■.托…t=0,采用R-F将非正念变量当罱正态化,得到等效正态分布的均值o:和标准差虹及可靠度指标B,由可靠度指标B的几何意义知。o;辟
开始时验算点未知,把6看成极限状态曲面上点P(■,爿:---37,的函数,通过优化求解,找到B最小值。求解可靠皮指标aJ以归结为以下约束优化模型: rain睁喜t华,2 s.,.Z=g(工i,x2’,…,工:=0 如极限状态方栉巾某个变最(X。可用其他变量表示,则上述模型jfIJ‘转化为无约束优化模型: 。。B!:手f生丛r+阻:坚:坠:盐尘}二剐 t∞oY?’【叫,J 3用MATLAB实现结构可靠度计算 3.1Matlab简介 Matlab是++种功能强、效率高、便.丁.进行科学和工程计算的交互式软件包,汇集了人量数学、统计、科学和工程所需的函数,MATI.AB具有编程简甲直观、用户界mf友善、开放性强等特点。将MATLAB用于蒙特卡罗法的一个显著优点是它拥有功能强大的随机数发生器指令。 3.2算例 3.2.I例:已知非线形极限状态方程z=g(t r'H=567f r-0.5H2=0’f、r服从正态分布。IIf=0.6,o r=0.0786;la|_ 2.18,o r_0.0654;H服从对数正态分布。u H= 3218,O。 =0.984。f、r、H相互独立,求可靠度指标B及验算点(,,r’,H‘。 解:先将H当量正念化:h=ln H服从正态分布,且 ,‘-““了:等专虿’=,。49?口二-、『五ir面_。。3
《计算方法》实验报告 指导教师： 学院： 班级： 团队成员：
一、题目 例2.7应用Newton 迭代法求方程210x x --=在1x =附近的数值解 k x ，并使其满足8110k k x x ---< 原理： 在方程()0f x =解的隔离区间[],a b 上选取合适的迭代初值0x ，过曲线()y f x =的点()() 00x f x ，引切线 ()()()1000:'l y f x f x x x =+- 其与x 轴相交于点：()() 0100 'f x x x f x =-，进一步,过曲线()y f x =的 点()()11x f x , 引切线 ()()()2111: 'l y f x f x x x =+- 其与x 轴相交于点：() () 1211 'f x x x f x =- 如此循环往复，可得一列逼近方程()0f x =精确解*x 的点 01k x x x ,,,,,其一般表达式为： ()() 111 'k k k k f x x x f x ---=- 该公式所表述的求解方法称为Newton 迭代法或切线法。
程序： function y=f(x)%定义原函数 y=x^3-x-1; end function y1=f1(x0)%求导函数在x0点的值 syms x; t=diff(f(x),x); y1=subs(t,x,x0); end function newton_iteration(x0,tol)%输入初始迭代点x0及精度tol x1=x0-f(x0)/f1(x0);k=1;%调用f函数和f1函数 while abs(x1-x0)>=tol x0=x1;x1=x0-f(x0)/f1(x0);k=k+1; end fprintf('满足精度要求的数值为x(%d)=%1.16g\n',k,x1); fprintf('迭代次数为k=%d\n',k); end 结果：
作者：夏云木子 1、 >> syms re(x) re(y) re(z) >> input('计算相对误差：'),re(x)=10/1991,re(y)=0.0001/1.991,re(y)=0.0000001/0.0001991 所以可知re(y)最小，即y精度最高 2、 >> format short,A=sqrt(2) >> format short e,B=sqrt(2) >> format short g,C=sqrt(2)
>> format long,D=sqrt(2) >> format long e,E=sqrt(2) >> format long g,F=sqrt(2) >> format bank,H=sqrt(2) >> format hex,I=sqrt(2) >> format +,J=sqrt(2) >> format,K=sqrt(2)
3、 >> syms A >> A=[sqrt(3) exp(7);sin(5) log(4)];vpa(pi*A,6) 4、1/6251-1/6252=1/6251*6252 5、（1）1/（1+3x）-(1-x)/(1+x)=x*(3*x-1)/[(1+3*x)*(1+x)] (2) sqrt(x+1/x)-sqrt(x-1/x)=2/x/[sqrt(x-1/x)+sqrt(x+1/x)] (3) log10(x1)-log(x2)=log10(x1/x2) (4) [1-cos(2*x)]/x =x^2/factorial(2)-x^4/factorial(4)+x^6/factorial(6)-…
function [x0,k]=bisect1(fun1,a,b,ep) if nargin<4 ep=1e-5; end fa=feval(fun1,a); fb=feval(fun1,b); if fa*fb>0 x0=[fa,fb]; k=0; return; end k=1; while abs(b-a)/2>ep x=(a+b)/2; fx=feval(fun1,x); if fx*fa<0 b=x; fb=fx; else a=x; fa=fx;
end end x0=(a+b)/2; >> fun1=inline('x^3-x-1'); >> [x0,k]=bisect1(fun1,1.3,1.4,1e-4) x0 = 1.3247 k = 7 >> 简单迭代法 function [x0,k]=iterate1(fun1,x0,ep,N) if nargin<4 N=500; end if nargin<3 ep=1e-5; end x=x0; x0=x+2*ep;
while abs(x-x0)>ep & k