MATLAB计算机仿真上机作业.pdf
杭电 matlab 上机实验
第一次上机内容：习题2：A=cos(pi/3)-(9-sqrt(2))^(1/3)习题11:x=0:2*pi/124:2*pi;y=cos(x)*(0.5+3*sin(x)/(1+x.*x));plot(x,y);习题13：x=-2:0.1:2;y=-2:0.1:2;[X,Y]=meshgrid(x,y);Z=X.*X.*exp(-X.^2-Y.^2);mesh(Z);习题15:N=2;t=0:0.1:2*pi;x=cos(t);y=sin(N*t+0);figure(1);subplot(2,2,1);plot(x,y);title('a=0');y=sin(N*t+pi/3);subplot(2,2,2);plot(x,y);title('a=pi/3');y=sin(N*t+pi/2);subplot(2,2,3);plot(x,y);title('a=pi/2');y=sin(N*t+pi);subplot(2,2,4);plot(x,y);title('a=pi');3.选择一幅彩色图片,将其变作灰度图，再做上下颠倒、左右颠倒、逆时针90，顺时针旋转180度A=imread('1.jpg','jpg');A1=double(A);A1=[A1(:,:,1)+A1(:,:,2)+A1(:,:,3)]/3;imshow(uint8(A1));[h,w]=size(A1);B=A1;for i=1:hB(i,:)=A1(h+1-i,:);endimshow(uint8(B));C=A1;for j=1:wC(:,j)=A1(:,w+1-j);endimshow(uint8(C));D=eye(w,h);for j=1:wD(j,:)=A1(:,w+1-j);endimshow(uint8(D));E=A1;for j=1:wfor i=1:hE(i,j)=A1(h+1-i,w+1-j);endendimshow(uint8(E));第二次上机内容1.鼠标点击空白图面，实时显示鼠标点和轨迹figure(1);axis([-1,1,-1,1]);[x,y,b]=ginput(1);plot(x,y,'+r');string=sprintf('(%5.3f,%5.3f)',x,y); text(x,y,string);hold on;while(1)if b==3break;endx1=x;y1=y;[x,y,b]=ginput(1);plot(x,y,'+');string=sprintf('(%5.3f,%5.3f)',x,y);text(x,y,string);plot([x1 x],[y1 y],'+r-');end2.正弦的轨迹动态表示axis([0,4*pi,-1,1]);x=0;y=sin(x);hold on;x1=x;y1=y;for x=pi/24:pi/24:4*pipause(0.1);y=sin(x);plot([x1 x],[y1 y]);x1=x;y1=y;end第三次上机内容：制作一界面，上有二个输入框、一个计算按钮、一个显示框、一个下拉选择框，做二数的可选择的四则运算。
计算机仿真技术-上机实验答案
计算机仿真技术-上机实验1. 已知传递函数为008.51026.2503247.5008.504.25)()()(23++++==s s s s s U s Y s G在MATLAB 中将传递函数转化为零极点和状态方程两种模型形式，要求写出MATLAB 指令和得到的模型方程.MATLAB 指令为：num=[25.04 5.008] ;den=[1 5.03247 25.1026 5.008] ;sys_tf=tf(num,den);[A,B,C,D]=tf2ss(num,den);sys_ss=ss(A,B,C,D)[z,p,k]=tf2zp(num,den);sys_zpk=zpk(z,p,k)得到：零极点模型为Zero/pole/gain:25.04 (s+0.2)---------------------------------(s+0.2078) (s^2 + 4.825s + 24.1)状态空间模型为：a =x1 x2 x3x1 -5.032 -25.1 -5.008x2 1 0 0x3 0 1 0b =u1x1 1x2 0x3 0c =x1 x2 x3y1 0 25.04 5.008d =u1y1 0Continuous-time model.2. 已知状态空间模型为2u-2y 形式：⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡--=⎥⎦⎤⎢⎣⎡21212121211100001001101142510u u x x y y u u x x x x在MATLAB 中将状态方程转化为传递函数，要求写出MATLAB 指令和得到的模型方程.注意：2个输入和2个输出，得到4个传递函数：G11=y1/u1, G12=y1/u2,G21=y2/u1,G22=y2/u2.MATLAB 指令为：A=[0 1 ; -25 -4] ; B=[1 1 ; 0 1] ; C=[1 0 ; 0 1] ; D=[0 0 ; 0 0] ;sys_ss=ss(A,B,C,D) ;[num,den]=ss2tf(A,B,C,D,1);sys_tf11=tf(num(1,:),den)sys_tf21=tf(num(2,:),den)[num,den]=ss2tf(A,B,C,D,2);sys_tf12=tf(num(1,:),den)sys_tf22=tf(num(2,:),den)Transfer function y1/u1s + 4--------------s^2 + 4 s + 25Transfer function: y2/u1-25--------------s^2 + 4 s + 25Transfer function: y1/u2s + 5--------------s^2 + 4 s + 25Transfer function: y2/u2s - 25--------------s^2 + 4 s + 253. 已知连续时间传递函数为62551010)()()(23++++==s s s s s U s Y s G 取采样时间为0.1s ，求离散化的传递函数，要求写出MATLAB 指令和离散化传递函数。
Matlab上机作业部分参考答案.ppt
[rank(A), rank([A B])]
ans =
34 由得出的结果看，A, [A;B] 两个矩阵的秩不同，故方程是
矛盾方程，没有解。
5. 试求下面齐次方程的基础解系
7. 建立如下一个元胞数组，现在要求计算第一个元胞第4行第 2列加上第二个元胞+第三个元胞里的第二个元素+最后一个元 胞的第二个元素。
a={pascal(4),'hello';17.3500,7:2:100}
解： >> a={pascal(4),'hello';17.3500,7:2:100} a=
[ 173/34, 151/34]
6. 求解方程组的通解
x1 2x2 4x3 6x4 3x5 2x6 4 2x1 4x2 4x3 5x4 x5 5x6 3
3x1 6x2 2x3 5x5 9x6 1 2x1 3x2 4x4 x6 8
4x2
5x3
2x4
x5
参考答案： （1） >> limit(sym('(tan(x) - sin(x))/(1cos(2*x))')) ans = 0 （2） >> y = sym('x^3 - 2*x^2 + sin(x)'); >> diff(y) ans = 3*x^2-4*x+cos(x) （3） >> f = x*y*log(x+y); >> fx = diff(f,x) fx = y*log(x+y)+x*y/(x+y)
MATLAB上机内容及作业
MATLAB上机内容及作业无约束优化求解函数fminsearch和fminunc求解无约束非线性优化问题的函数有fminsearch 函数、fminunc 函数。
函数fminsearch和fminunc功能相同，但fminunc函数可以得到目标函数在最优解处的梯度和Hessian矩阵值。
无约束优化数学模型为：min f(X) X∈R n求解无约束非线性优化问题的步骤为：第一步：先编写目标函数的M文件；第二步：再在命令窗口中调用相应的优化函数。
1、fminsearch函数调用格式为[x, fval]=fminsearch(@myfun, x0)输出参数的含义：x：返回最优解的设计变量的值；fval：在最优设计变量值时，目标函数的最小值；exitflag：返回算法终止的标志，有以下几种情况，＞0 表示算法收敛于最优解处；＝0 表示算法已经达到迭代的最大次数；＜0 表示算法不收敛。
output：返回优化算法信息的一个数据结构，有以下信息：output.iteration 表示迭代次数output.algorithm 表示所采用的算法output.funcCount 表示函数评价次数输入参数的含义：@myfun：目标函数的M文件，在其前要加“@”，或表示为'myfun' ，myfun自己可以任意命名；x0：在调用该优化函数时，需要先对设计变量赋一个初始值；2、fminunc函数的调用格式[x, fval]=fminunc (@myfun, x0)grad：返回目标函数在最优解处的梯度信息；hessian：返回目标函数在最优解处的hessian矩阵信息。
其余含义同上。
3、实例已知某一优化问题的数学模型为：min f(X)=3x12+2x1x2+x22X∈R n用MA TLAB程序编写的代码为：第一步：首先编写目标函数的.m文件，并保存为examplefsearch.m文件（先单击file菜单，后点击New 命令中的M—file，即可打开M文件编辑窗口进行代码的编辑，在英文状态下输入程序代码），代码为：function f=examplefsearch（x）f=3*x(1)^ 2+2*x(1)*x(2)+x(2)^2;第二步：在Command窗口中调用fminsearch函数，代码为：x0=[1;1]; %赋初值[x,fval]=fminsearch(@examplefsearch,x0) %回车即可调用fminsearch函数，得到结果输出最优解结果为：x=1.0e-0.08* -0.7914 0.2260 %分别为x1和x2的最优点的值（近似为0）fval=1.5722e-016 %对应最优点的最优目标函数值（近似为0）4、作业已知几个优化问题的数学模型分别为：（1）min f(X)=0.1935x1 x22 x32(4+6x4) X∈R4（2）min f(X)= (x13+cos x2+log x3）/ e x1 X∈R3（3）min f(X)=2x13+4x1x23 -10x1x2+x33X∈R3试用MATLAB编程分别求解上述优化问题的最优解。
计算机仿真上机
实验一：MATLAB基础入门实验目的：熟悉MATLAB环境，掌握一维数组的创建，二维数组的创建。
(1)a=[1,2,pi,9,0]；b=0:2:10；c=linspace(1,2,10)a(2) b(5) c(6)(2)二维数组的创建：a=[1,2,3,0;9,22,1,1];a(1,2) a(2,3) a(:) a(:,:) a(:,1) a(2,:)[1,2;3,4]+10；[1,2;3,4]*[0.1,0.2;0.3,0.4]；[1,2;3,4].\[20,10;9,2]sin([1,2;3,4])[a,b]=find([1,2;3,4]>=[4,3;2,1]) [a,b]find([1,2;3,4]>=[4,3;2,1])(3) 观察下列计算结果，理解这些命令的意义Clear；Which ；clc；的作用实验二：MATLAB基础入门实验目的：掌握MATLAB中基本的二维绘图t=0:pi/20:2*pi;y=sin(t);plot(t,y)t=[0 1]; x=[1 2]; y=[x;3 4]; z=[y;5 6];plot(t,x,'r') plot(t,y, 'b') plot(t,y') plot(t,z) plot(t,z')线型-实线:虚线-.点划线--间断线点标记.点o小圆圈x叉子符+加号*星号s方格d菱形^朝上三角v朝下三>朝右三角<朝左三角p五角星h六角星颜色y黄色m品红色c青色r红色g绿色b蓝色w白色上机指出以下各个绘图命令的输出图形分别是什么t=0:pi/20:2*pi; plot(t,sin(t) ,'r:>');hold on;plot(t,cos(t),'b-.h')title('sin(t),cos(t)的函数图形');xlabel('t=0:pi/20:2*pi;')ylabel('sin(t),cos(t)');legend('sin(t)','cos(t)')体验grid on /grid off;hold on/hold off;figure(2)的作用利用plot 函数在一个坐标系下绘制以下函数的图形：y1=sin(x),y2=cos(x),y3=sin(2*x),x在0到2*pi区间y1 用黑色间断线点标记为星号y2 用红色实线点标记为小方格y3 用蓝色虚线点标记为小圆圈实验三：MATLAB编程入门实验目的：掌握matlab编程的基本知识。
Matlab上机实验答案.docx
实验一MATLAB 运算基础Matlab 上机实验答案1. 先求下列表达式的值，然后显示MATLAB X 作空间的使用情况并保存全部变量。
» z1 =2*sin(85*pi/180)/(1 +exp(2)) z1 =0.2375⑵ z? =*ln(x +Jl + x?),其中 「45 » x=[2 1+2i;・0.45 5];» z2=1/2*log(x+sqrt(1 +x A2)) z2 =0.7114 ・ 0.0253i 0.8968 + 0.3658i 0.2139+ 0.9343i1.1541 ・0.0044ie 0.3asin(a+0.3) + ln^^a = —3Q — 29 2.9,3.0(1)2sin85°» a=-3.0:0.1:3.0;» z3=(exp(0.3.*a)-exp(-0.3.*a))./2.*sin(a+0.3)+log((0.3+a)./2)(» z33=(exp(0.3*a)-exp(-0.3*a))/2.*sin(a+0.3)+log((0.3+a)/2)可以验证z3==z33,是否都为1)z3 =Columns 1 through 50.7388 + 3.1416i 0.7696 + 3.1416i 0.7871 + 3.1416i 0.7913 + 3.1416i 0.7822 + 3.1416iColumns 6 through 100.7602 + 3.1416i 0.7254 + 3.1416i 0.6784 + 3.1416i 0.6196 + 3.1416i 0.5496 + 3.1416iColumns 11 through 150.4688 + 3.1416i 0.3780 + 3.1416i 0.2775 + 3.1416iColumns 16 through 20-0.0771 + 3.1416i -0.2124 + 3.1416i -0.3566 + 3.1416i -0.5104 + 3.1416i -0.6752 + 3.1416iColumns 21 through 25-0.8536 + 3.1416i -1.0497 + 3.1416i -1.2701 + 3.1416i -1.5271+3.1416i -1.8436 + 3.1416iColumns 26 through 30-2.2727 + 3.1416i -2.9837 + 3.1416i -37.0245 + O.OOOOi-3.0017 + O.OOOOi -2.3085 + O.OOOOiColumns 31 through 35-1.8971 + O.OOOOi -1.5978 + O.OOOOi -1.3575 + O.OOOOi -1.1531 + O.OOOOi -0.9723 + O.OOOOiColumns 36 through 400.4841 + O.OOOOi0.6474 + O.OOOOi0.6119 + O.OOOOi 0.5777 + O.OOOOi 0.5327 + O.OOOOi-0.8083 + O.OOOOi -0.6567 + O.OOOOi -0.5151 + O.OOOOi-0.3819 + O.OOOOi -0.2561 + O.OOOOiColumns 41 through 45-0.1374 + O.OOOOi -0.0255 + O.OOOOi 0.0792 + O.OOOOi0.1766 + O.OOOOi 0.2663 + O.OOOOiColumns 46 through 500.3478 + O.OOOOi 0.4206 + O.OOOOi0.5379 + O.OOOOi 0.5815 + O.OOOOiColumns 51 through 550.6145 + O.OOOOi 0.6366 + O.OOOOi0.6470 + O.OOOOi 0.6351 + O.OOOOiColumns 56 through 600.4774 + O.OOOOi 0.4126 + O.OOOOiColumn 610.3388 + O.OOOOir 0<r<l(4) z4=< t2 -1 l<t<2 9其中Z=0:0・5:2・ 5尸—2r + l 2<r<3»t=0:0.5:2.5;»z4=(t>=0&t<1 ).*(t.A2)+(t>=1 &t<2).*(t.A2-1 )+(t>=2&t<3).*(t.A2-2.*t+1)z4 =0 0.2500 0 1.2500 1.0000 2.25002.已知：求下列表达式的值：(1) A+6*B和A・B+I (其中I为单位矩阵)(2) A*B 和 A.*B⑶A八3和A/3⑷A/B及B\A(5)[A,B]和[A([1,3],:);B^2]» A=[12 34 -4;34 7 87;3 65 7];» B=[1 3-1;2 0 3;3 -2 7];» A+6*Bans =» A.*Bans =12 102 468 0 2619 ・130 49» A A3 ans =37226 23382424737078688 45414214918848604600766118820»A.A3ans =» B\A ans =ans =12 3 4 11 2067.0000-134.0000 68.00003.设有矩阵A和B(1) 求它们的乘积C。
Matlab上机作业部分参考答案
1
1
( x s inx )x2d x
的命
(2 ,4 .9 ) 4. 假设已知一组数据，试用插值方法绘制出 x
区间内的光滑函数曲线，比较各种值算法的优劣。
x
-2
i
-1.7
-1.4
-1.1
-0.8
-0.5
-0.2
0.1
0.4
0.7
1
1.3
y
x
i
.102 .117 .131 .144 .156 .166 .173 .177 .178 .176 .171 .163 89 41 58 83 56 22 32 5 53 35 09 02
11. 用单选框做一个如图所示的界面，通过选择不同的单选 框来决定使用不同的色彩图。
1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8
0
2
4
6
8
10
12
14
假设 求
f (x ,y ) e d t
t2 0
x y
x 2 f 2 f 2 f 2 ，试 2 y x x y y2
.17853 .17635 .17109 .16302
y
x
i
1.6
i
1.9
2.2
2.5
2.8
3.1
3.4
3.7
4
4.3
4.6
4.9
.15255 .1402
.12655 .11219 .09768 .08353 .07015 .05786 .04687 .03729 .02914 .02236
y
i
【求解】用下面的语句可以立即得出给定样本点数据的三次插值与样条 插值，得出的结果如，可见，用两种插值方法对此例得出的结果几乎一 致，效果均很理想。 >> x=[-2,-1.7,-1.4,-1.1,-0.8,-0.5,-0.2,0.1,0.4,0.7,1,1.3,... 1.6,1.9,2.2,2.5,2.8,3.1,3.4,3.7,4,4.3,4.6,4.9]; y=[0.10289,0.11741,0.13158,0.14483,0.15656,0.16622,0.17332,... 0.1775,0.17853,0.17635,0.17109,0.16302,0.15255,0.1402,... 0.12655,0.11219,0.09768,0.08353,0.07019,0.05786,0.04687,... 0.03729,0.02914,0.02236]; x0=-2:0.02:4.9; y1=interp1(x,y,x0,'cubic'); y2=interp1(x,y,x0,'spline'); plot(x0,y1,':',x0,y2,x,y,'o')
matlab上机练习（附答案）
1.以下两种说法对吗？（1）MATLAB进行数值的表达精度与其指令窗中的数据显示精度相同。
（2）MATLAB指令窗中显示的数据有效位数不超过七位。
2.历史指令窗所记录的内容与diary指令所产生的“日志”内容有什么不同？DIARY filename causes a copy of all subsequent command window inputand most of the resulting command window output to be appended to thenamed file. If no file is specified, the file 'diary' is used.DIARY OFF suspends it.DIARY ON turns it back on.DIARY, by itself, toggles the diary state.Use the functional form of DIARY, such as DIARY('file'),when the file name is stored in a string.3.如何把用户自己的“工作目录”永久地设置在MATLAB的搜索路径上？“位于搜索路径上的目录”与“当前目录”在MATLAB中的功用相同吗？4.如何向MATLAB工作空间输入一个含有100个左右元素的一维或二维数值数组？用直接键入法？用数组编辑器？用M文件编辑器？5.运用数组算术运算符去掉下面程序里的for/end循环：x=11:15for k=1:length(x)z(k)=x(k)^2+2.3*x(k)^0.5;endx=11:15 z1=x.^2+2.3*x.^0.56.不使用数组算术运算符，重写下面的程序代码：x=[2 1 4]z=1./(1+x.^2)x=2;k=1;while i<=4,z2(k)=1/(1+i^2);i=i+1;x=x+1;end7．某公司销售电脑打印机的价格方案如下：（）如果顾客只买一台打印机，则一台的基本价格为$150。
matlab上机练习题答案(可编辑修改word版)
a ⎣ ⎦ 1. 计算 a = 5⎦ >> a=[6 9 3;2 7 5]; >> b=[2 4 1;4 6 8]; >> a.*b ans =⎣4 8⎦>> d=deconv([3 13 6 8],[1 4]) d =312⎡2 4 7 4⎤⎡8⎤1236 36 求欠定方程组⎢9 3 5 6⎥ x = ⎢5⎥ 的最小范数解84240⎡4 9 2⎤⎡37⎤⎣ >> a=[2 4 7 4;9 3 5 6]; >> b=[8 5]'; ⎦ ⎣ ⎦2. 对于 AX = B ，如果 A = ⎢7 6 4⎥ ， B = ⎢26⎥ ，求解 X 。
>> x=pinv(a)*b⎢ ⎥ ⎢ ⎥x =>> A=[4 9 2;7 6 4;3 5 7]; >> B=[37 26 28]’; >> X=A\B X =⎢⎣3 5 7⎥⎦ ⎢⎣28⎥⎦-0.2151 0.4459 0.7949 0.27077 用符号函数法求解方程 a t 2+b*t +c=0-0.5118 >> r=solve('a*t^2+b*t+c=0','t') 4.0427 r =1.3318[ 1/2/a*(-b+(b^2-4*a*c)^(1/2))] ⎡1 2 5 ⎤ ⎡8 - 7 4⎤ [ 1/2/a*(-b-(b^2-4*a*c)^(1/2))]3. a = ⎢3 6 - 4⎥ ， b = ⎢3 6 2⎥ ，观察 a 与 b 之间的⎣ ⎦ ⎣ ⎦⎡a 11 a 12 ⎤六种关系运算的结果 >> a=[1 2 3;4 5 6];8 求矩阵 A = ⎢ ⎣ 21 ⎥ 的行列式值、逆和特征根a 22 ⎦>> b=[8 –7 4;3 6 2];>> syms a11 a12 a21 a22;>> a>b >> A=[a11,a12;a21,a22] ans =>> AD=det(A) % 行列式 0 1 0 >> AI=inv(A) % 逆 11>> AE=eig(A) % 特征值 >> a>=b ans =0 1 0 1 01>> a<b ans =1 0 1 0 1>> a<=b ans =1 0 1 010 A = [ a11, a12][ a21, a22] AD =a11*a22-a12*a21 AI =[ -a22/(-a11*a22+a12*a21), a12/(-a11*a22+a12*a21)] [ a21/(-a11*a22+a12*a21), -a11/(-a11*a22+a12*a21)] AE = [1/2*a11+1/2*a22+1/2*(a11^2- 2*a11*a22+a22^2+4*a12*a21)^(1/2)] [1/2*a11+1/2*a22-1/2*(a11^2->> a==b2*a11*a22+a22^2+4*a12*a21)^(1/2)] ans =9 因式分解： x 4 - 5x 3 + 5x 2 + 5x - 60 0 0 >> syms x;0 00 >> f=x^4-5*x^3+5*x^2+5*x-6; >> a~=b ans =>> factor(f) ans =(x-1)*(x-2)*(x-3)*(x+1)4 计算多项式乘法(x 2+2x +2)(x 2+5x +4) ⎡10 f = ⎢ x 21 x ⎤⎥ ，用符号微分求 df/dx 。
Matlab上机操作作业(4)参考答案
上机操作（4）内容： M 文件编程一、结合教材第6章内容以及实验指导书内容掌握一般M 脚本和函数文件的编写。
二、 完成以下习题，并按照要求将作业提交。
（1） 将调试通过的程序序列复制到作业题目下方，并将文件保存为word 文档，文件名为：姓名+班级学号末两位数 (如：王丹妮9301)（2） 作业提交方式：A 现场网上提交（鼠标右点击电脑右下方的人头或燕子标记，按照菜单操作即可。
B 若不能提交，发送至我的邮箱363586609@ .发送时邮件标题为：文件名作业4（如：王丹妮9301作业4）1. 编写M 脚本文件，画出下列分段函数所表示的曲面。
（30分）⎪⎪⎩⎪⎪⎨⎧-≤+≤+<->+=+-------15457.0117575.015457.0),(215.175.375.0216215.175.375.02112122212212122x x e x x e x x e x x p x x x x x x x x %参考答案a=2;b=2;x=-a:0.1:a;y=-b:0.1:b;for i=1:length(y)for j=1:length(x)if x(j)+y(i)>1z(i,j)=0.5457*exp(-0.75*y(i)^2-3.75*x(j)^2-1.5*x(j));elseif x(j)+y(i)<=-1z(i,j)=0.5457*exp(-0.75*y(i)^2-3.75*x(j)^2+1.5*x(j));else z(i,j)=0.7575*exp(-y(i)^2-6.*x(j)^2);endendendaxis([-a,a,-b,b,min(min(z)),max(max(z))]);surf(x,y,z);xlabel('x1');ylabel('x2');zlabel('p(x1,x2)');2．编写一个M 函数文件，该文件应有函数声明行、帮助说明和程序编写人姓名。
Matlab第二、三次上机作业-推荐下载
第二次上机作业一. 任务： 用MATLAB 语言编写连续函数最佳平方逼近的算法程序（函数式M 文件）。
并用此程序进行数值试验，写出计算实习报告。
二. 程序功能要求：1.用Lengendre 多项式做基，并适合于构造任意次数的最佳平方逼近多项式。
可利用递推关系0112()1,()()(21)()(1)()/2,3,.....n n n P x P x x P x n xP x n P x nn --===---⎡⎤⎣⎦=2.程序输入：（1）待求的被逼近函数值的数据点(可以是一个数值或向量)0x （2）区间端点：a,b 。
3. 程序输出：（1）拟合系数：012,,,...,nc c c c （2）待求的被逼近函数值00001102200()()()()()n n s x c P x c P x c P x c P x =++++ 三：数值试验要求：1.试验函数：；也可自选其它的试验函数。
()cos ,[0,4]f x x x x =∈+2.用所编程序直接进行计算，检测程序的正确性，并理解算法。
3.分别求二次、三次、。
最佳平方逼近函数。
）x s (4.分别作出逼近函数和被逼近函数的曲线图进行比较。
）x s ()(x f （分别用绘图函数plot(,s())和fplot(‘x cos x ’,[x 1 x 2,y 1,y 2])）0x 0x 解题思路：参照应用数值分析书P259“算法7-1”，利用Legendre 多项式对f(x)∈C(a,b)的最佳平方逼近写出以下算法：M 文件1：%文件名：GLAppro.mfunction[poly,yy,delta]=GLAppro(f,n,a,b,xx)%功能：利用Gauss Legendre 多项式求函数的最佳平方逼近%输入：f ——被逼近函数；a,b ——逼近区间；xx ——欲求的逼近点%n ——逼近的L 多项式的次数（标量时为最高次数，向量时为其所选择的的逼近次数）%输出：poly ——所求的逼近多项式系数（降序）；yy ——逼近店的值；delta ——逼近误差N=numel(n);if N>1id=n+1;elseN=n+1;enddelta=quad(@myfun,-1,1,[],[],f,a,b);c=zeros(1,N);poly=zeros(1,id(N)); for k=1:Nc(k)=(2*id(k)-1)*quad(@fLegPoly,-1,1,[],[],f,id(k)-1,a,b)/2;delta=delta-c(k)^2*2/(2*id(k)-1);endif nargin==5t0=(2*xx-a-b)/(b-a);yy=zeros(size(xx));for k=1:Np=LegPoly(id(k)-1);yy=yy+c(k)*polyval(p,t0);poly(id(N)-id(k)+1:(id(N)))=poly(id(N)-id(k)+1:(id(N)))+c(k)*p;endelsefor k=1:Np=LegPoly(k-1);poly(N-k+1:N)=poly(N-k+1:N)+c(k)*p;endendM文件2：function y=myfun(t,f,a,b)%功能：GLAppro子函数，变换到区间[-1,1]x=(b-a)*t/2+(b+a)/2;y=f(x).*f(x);M文件3：function f=fLegPoly(t,f1,n,a,b)%功能：GLAppro子函数，求变换后的积分函数p=LegPoly(n);x=(b-a)*t/2+(b+a)/2;f=f1(x).*polyval(p,t);M文件4：function p=LegPoly(n)%功能：递归法求n次Gauss Legendre多项式%输入：n——多项式次数%输出：p——降幂排列的多项式系数p0=1;p1=[1,0];if n==0p=p0;elseif n==1p=p1;elsep=((2*n-1)*[LegPoly(n-1),0]-(n-1)*[0,0,LegPoly(n-2)])/n;end%文件名：homework2.mf=@(x)(x.*cos(x));n=cell(3,1);n{1}=2;n{2}=3;n{3}=4;a=0;b=4;x0=linspace(a,b)';color=['k','g','b'];y=f(x0);plot(x0,y,'r-','linewidth',1.5);hold on;syms t x;for i=1:3[poly,py,delta]=GLAppro(f,n{i},a,b,x0);pt=vpa(poly2sym(poly,t),4);poly=simple(subs(pt,t,(2*x-a-b)/(b-a)));poly=vpa(poly,4);disp(['所求的最佳评分那个逼近多项式为（n=' int2str(n{i}) '):']);disp(poly);disp(['误差为:' num2str(delta)]);plot(x0,py,color(i),'linewidth',1.5)end分别求两次、三次和四次最佳平方逼近函数s(x)输入：homework2输出：所求的最佳评分那个逼近多项式为（n=2):- 0.1599*x^2 - 0.572*x + 0.8268误差为:0.45032所求的最佳评分那个逼近多项式为（n=3):0.3818*x^3 - 2.45*x^2 + 3.093*x - 0.3948误差为:0.02393所求的最佳评分那个逼近多项式为（n=4):0.08539*x^4 - 0.3014*x^3 - 0.6939*x^2 + 1.531*x - 0.08255误差为:0.002258调整图象并标注：yx第三次上机作业一. 任务一：用MATLAB 语言编写化工原理流体力学中d-u-Re-λ试差法的算法程序（函数式M 文件）。
计算机仿真(Matlab)上机报告
计算机仿真上机报告班级：信息101学号：201027012 姓名：张化迪指导老师：李贺实验一1、实验要求：安装Matlab软件，熟悉Matlab语言的基本语法格式及用法2.实验过程：3、实验感受：本次实验内容主要是安装Matlab以及熟悉Matlab的界面和基本操作，安装Matlab软件后，按照老师的要求执行了Matlab命令。
我熟悉了Matlab软件的安装步骤并了解了Matlab的基本语句和用法。
实验二1、实验要求：用Matlab编程语言实现AMI编码，要求可以将文本文件中的数据读入matlab程序，然后转换为数值量进行AMI编码，然后存入一个新的文本文件中。
并且将编码前后的数据用画图函数在程序中实现，运行并测试。
2、实验程序：3、实验结果4.实验分析本次实验中要求实现AMI的编码转码和图像表示，在本次实验中，在工作区中建立对应的AMI的编码的文件，注意码元1 、0 之间要有空格，这样可以确保Matlab程序对文件的读入和转换是正确的。
通过对转码后的波形可得出，利用上述程序成功实现了对AMI码的转换,同时另一个生成的文本文档也验证了该程序实现老师要求的AMI编码的可行性。
实验三1、实验要求：以经典的二阶电路为例（例如RLC电路），利用数值分析函数（例如龙格库塔函数）用Matlab语言编写程序实现其参数的计算及观测2、程序：方程定义的M文件（文件1）使用ode32仿真的M文件（文件2）t0=0;3、仿真结果：4、实验总结：本实验通过Matlab运用数值分析函数对二阶电路进行分析，对二阶电路各参数进行计算，得到相应图像。
通过此次实验我也掌握了Matlab中M格式文件的用法，并熟悉了数值分析函数的使用方法以及使用Matalab对二阶电路编程分析的方法。
实验四1、实验要求：熟悉Simulink建模及通信、数字信号处理工具箱，建立简单的数字信号BPSK 的调制解调模型，并借助示波器等虚拟工具观测数据的正确性。
matlab上机考试题及答案
matlab上机考试题及答案1. 题目：编写一个MATLAB函数，计算并返回一个向量中所有元素的平方和。
答案：函数定义如下：```matlabfunction sumOfSquares = calculateSumOfSquares(vector)sumOfSquares = sum(vector.^2);end```2. 题目：使用MATLAB的内置函数，找出一个矩阵中的最大元素及其位置。
答案：可以使用`max`函数来找出矩阵中的最大元素，同时使用`find`函数来获取其位置。
示例代码如下：```matlabA = [1, 2, 3; 4, 5, 6; 7, 8, 9];[maxValue, linearIndex] = max(A(:));[row, col] = ind2sub(size(A), linearIndex);```3. 题目：给定一个向量，使用MATLAB编写代码，实现向量元素的逆序排列。
答案：可以使用`flip`函数来实现向量的逆序排列。
示例代码如下：```matlabvector = [1, 2, 3, 4, 5];reversedVector = flip(vector);```4. 题目：编写一个MATLAB脚本，计算并绘制一个正弦波的图像。
答案：可以使用`sin`函数生成正弦波数据，并使用`plot`函数绘制图像。
示例代码如下：```matlabx = linspace(0, 2*pi, 100);y = sin(x);plot(x, y);xlabel('x');ylabel('sin(x)');title('Sine Wave');```5. 题目：给定一个3x3的矩阵，使用MATLAB编写代码，计算其行列式。
答案：可以使用`det`函数来计算矩阵的行列式。
示例代码如下：```matlabmatrix = [1, 2, 3; 4, 5, 6; 7, 8, 9];determinant = det(matrix);```结束语：以上是MATLAB上机考试的题目及答案，希望能够帮助大家更好地掌握MATLAB的编程技巧和函数使用。
《MATLAB及应用》上机作业-2015.doc
《MA TLAB及应用》上机作业学院名称：（四号宋体）专业班级：（四号宋体）学生姓名：（四号宋体）学生学号：（四号宋体）年月《MATLAB及应用》上机作业要求及规范一、作业提交方式：word文档打印后提交。
二、作业要求：1．封面：按要求填写学院、班级、姓名、学号，不要改变封面原有字体及大小。
2．内容：只需解答过程（结果为图形输出的可加上图形输出结果），不需原题目；为便于批阅，题与题之间应空出一行；每题答案只需直接将调试正确后的M文件内容复制到word 中（不要更改字体及大小），如下所示：%作业1_1clcA=[1 2 3 4;2 3 5 7;1 3 5 7;3 2 3 9;1 8 9 4];B=[1+4*i 4 3 6 7 7;2 3 3 5 5 4+2*i;2 6+7*i 5 3 4 2;1 8 9 5 4 3];C=A*BD=C(4:5,4:6)三、大作业评分标准：1．提交的打印文档是否符合要求；2．作业题的解答过程是否完整和正确；3．答辩过程中阐述是否清楚，问题是否回答正确；4．作业应独立完成，严禁直接拷贝别人的电子文档，发现雷同者都以无成绩论处。
作业11、用MATLAB 可以识别的格式输入下面两个矩阵12342357135732391894A ⎛⎫⎪⎪ ⎪= ⎪⎪ ⎪⎝⎭，144367723355422675342189543i i B i +⎛⎫⎪+⎪= ⎪+ ⎪⎪⎝⎭再求出它们的乘积矩阵C ，并将C 矩阵的右下角23⨯子矩阵赋给D 矩阵。
赋值完成后，调用相应的命令查看MATLAB 工作空间的占有情况。
2、设矩阵16231351110897612414152⎛⎫⎪⎪ ⎪ ⎪⎝⎭，求A ，1A -，3A ，12A A -+，1'3A A --，并求矩阵A 的特征值和特征向量。
3、解下列矩阵方程:010100143100001201001010120X -⎛⎫⎛⎫⎛⎫⎪ ⎪ ⎪=- ⎪ ⎪ ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ 4、求多项式322()(53513)(453)(313)f x x x x x x x =++++++当3x =时的值，并求()f x 的导数。
Matlab上机作业部分参考答案
3. 有一组学生的考试成绩如下表，根据规定，成绩在100分时 为满分，成绩在90~99之间时为优秀，成绩在80~89分之间时 为良好，成绩在60~79分之间为及格，成绩在60分以下时为不 及格，编制一个根据成绩划分等级的程序。
学生 姓名 王 张 刘 李 陈 杨 于 黄 郭 赵成 绩Βιβλιοθήκη 728356
94
100
z1 x2 2y 2 , 2. 列出求下列空间曲面
（提示find 命令）
z2 2 x 3 y 交线的程序
clc clear x=[-5:0.5:5]; [X,Y]=meshgrid(x); z1=X.^2-2*Y.^2; z2=X.*2-Y.*3; xlabel('x') ylabel('y') zlabel('z')
surf(X,Y,z1) hold on surf(X,Y,z2) k=find(abs(z1-z2)<0.5); x1=X(k) y1=Y(k) z3=x1.^2-2*y1.^2 hold on
88
96
68
54
65
4. 编制写一个函数实现，输入一个3x4矩阵，求其中最大元素 和最小元素（禁止直接用Matlab命令） 5. 编程，可输入有10个实数的数组，调用函数f求他们的平均值 V，然后输出其中与V之差的绝对值为最小的那个数组元素。 6. 编制写一个函数实现，一个给定二维（3x3）数组转置。 （禁止直接用Matlab命令）
7.
编制写一个函数实现，输入16进制，输出10进制。（禁止直 接用Matlab命令）
8. 编程序实现功能：在一个元素值按升序存放的整型数组中插 入一个数，使得插入后的数组元素仍然有序。 9. 编制写一个函数实现，用“冒泡法” 排列输入的10个数字。
计算机仿真实验指导书（第一、二次上机）20151130
计算机仿真实验指导书（第⼀、⼆次上机）20151130计算机仿真实验指导书⽬录第⼀章MATLAB及其⼯作环境介绍 (1)1.1 MATLAB简介 (1)1.2 MATLAB的⼯作环境介绍 (1)1.3 MATLAB的基本管理命令 (4)第⼆章MATLAB的数值计算功能 (6)2.1 变量与赋值语句 (6)2.2 MATLAB矩阵 (6)2.3 MATLAB表达式 (11)2.4 MATLAB常⽤数学函数 (12)2.5 矩阵的基本运算 (13)2.6 矩阵的.点操作运算 (17)2.7 多项式及其运算 (18)第三章MATLAB程序设计⼊门 (20)3.1 M⽂件 (20)3.2 数据的输⼊输出 (22)3.3 全局变量和局部变量 (24)3.4 程序流程控制 (24)第四章MATLAB的符号运算功能 (29)4.1 建⽴符号对象 (29)4.2 符号算术运算 (30)4.3 符号微积分运算 (33)4.4 符号函数的可视化 (35)第五章MATLAB的可视化功能 (38)5.1 ⼆维图形 (38)5.2绘制三维图形 (43)5.3 特殊坐标图形 (45)5.4 图形句柄 (46)第六章控制系统的频域与时域分析 (52)6.1线性系统的描述 (52)6.2 模型之间的转换 (56)6.3 时域响应分析 (59)6.4 频域响应分析 (66)6.5 根轨迹分析 (70)第⼀章MATLAB及其⼯作环境介绍1.1 MATLAB简介MATLAB是matrix和laboratory前三个字母的缩写，意思是实验室矩阵。
MATLAB语⾔是⼀种⼴泛应⽤于⼯程计算及数值分析领域的新型⾼级语⾔，⾃1984年由美国MathWorks公司推向市场以来，经过⼗多年的发展与完善，MATLAB已发展成为由MATLAB语⾔、MATLAB⼯作环境、MA TLAB图象处理系统、MATLAB 数学函数库和MATLAB应⽤程序接⼝五⼤部分组成的集数值计算、图形处理、程序开发为⼀体的功能强⼤的体系。
完整word版MATLAB上机作业
MATLAB 上机作业1对以下问题编写M文件：(1)用起泡法对10个数由小到大排序。
即将相邻两个数比较，将小的调到前头。
fun cti on f=qip aofa(x)for j=9:-1:1for i=1:jif(x(i)>x(i+1))t=x(i);x(i)=x(i+1);x(i+1)=t;end end end f=xx=rou nd(10*ra nd(1,10))qip aofa(x);(2)有一个4X 5矩阵，编程求出其最大值及其所处的位置。
fun cti on f=zuidazhi(x)a=1;b=1;c=x(1,1);for i=1:4for j=1:5 if x(i,j)>c a=i;b=j;c=x(i,j);endend end f=[c,a,b] x=ra nd(4,5) zuidazhi(x)20(3)编程求送n! ond：fun cti on f=qiuhe(x) b=0;for i=1:xa=p rod(1:i); b=b+a;end f=bqiuhe(20)(4)一球从100米高度自由落下，每次落地后反跳回原高度的一半，再落下。
求它在第10次落地时，共经过多少米？第10次反弹有多高？fun cti on f=gao(x)b=x;for i=2:10x=x/2;a=x*2; b=b+a;end f=[b x/2]gao(100)⑸有一函数f（x,y）=x2 +sin xy+2y，写一程序，输入自变量的值，输出函数值。
Function f=fun(x)f=x(1)A2+si n(x(1)*x(2))+2*x (2)MATLAB 上机作业21. 求和 Y =1 +4+42+440。
syms k s=4^k;S=symsum(s,k,0,40) 2. 求函数f （X ）=2x ' -6x 2-18x + 7的极值，并作图。
y='2*x^3-6*x^2-18*x+7'; y_='-2*xA3+6*xA2+18*x-7';[x_mi n, y_mi n]=fmi nbn d(y,-7,7) [x_max,y_max]=fmi nbn d(y_,-7,7) ezpl ot(y) 3. 设 y =e x sinx-7cosx+5x 2，求 业=?,3^ = dx dxy=ex p(x)*si n(x)-7*cos(x)+5*xA2; dy=diff(y)d2y=diff(y,2) 兀/2 ■ ------------4.求积分 h =J 0 J1-2si n2xdx 。
《MATLAB及其应用》上机作业.doc(修订于2009.11.19)
《MATLAB及其应用》上机作业学院名称：（四号宋体）专业班级：（四号宋体）学生姓名：（四号宋体）学生学号：（四号宋体）年月作业11.用MATLAB 可以识别的格式输入下面两个矩阵12342357135732391894A ⎛⎫⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭,144367723355422675342189543i i B i +⎛⎫ ⎪+⎪= ⎪+ ⎪⎪⎝⎭再求出它们的乘积矩阵C ,并将C 矩阵的右下角23⨯子矩阵赋给D 矩阵。
赋值完成后,调用相应的命令查看MATLAB 工作空间的占有情况。
2.设矩阵16231351110897612414152⎛⎫⎪⎪ ⎪ ⎪⎝⎭,求A ,1A -,3A ,12A A -+,1'3A A --,并求矩阵A 的特征值和特征向量。
3.解下列矩阵方程:010100143100001201001010120X -⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪=- ⎪ ⎪ ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭4.一球从100米高度自由落下,每次落地后反跳回原高度的一半,再落下。
求它在第10次落地时,共经过多少米?第10次反弹有多高?5.用MATLAB 语言实现下面的分段函数5,1(),25,y f x x ⎧⎪⎪==⎨⎪-⎪⎩101010x x x >≤<-6.分别用for 和while 循环编写程序,求出6323626312122222i i K ===++++++∑并考虑一种避免循环的简洁方法来进行求和,并比较各种算法的运行时间。
7.应用MA TLAB 语言及二分法编写求解一元方程32()1459700f x x x x =-+-=在区间[3，6]的实数解的算法，要求绝对误差不超过0.001。
8.二阶系统的单位阶跃响应为1cos )t y a ζζ-=+，在同一平面绘制ζ分别为0，0.3，0.5，0.707的单位阶跃响应曲线。
要求：(1)四条曲线的颜色分别为蓝、绿、红、黄，线型分别为“——”、“……”、“oooooo”、“++++++”； (2)添加横坐标轴和纵坐标轴名分别为“时间t”和“响应y”，并在平面图上添加标题“二阶系统曲线”和网格；(3)在右上角添加图例（即用对应的字符串区分图形上的线），并分别在对应的响应曲线的第一个峰值处标示“zeta ＝0”、“zeta ＝0.3”、“zeta ＝0.5”、“zeta ＝0.707”。
MATLAB计算机仿真上机作业
1、下列表达式表达了在主接触应力F 的作用下，两个球体被挤压到一起时在x,y,z 方向上所产生的接触应力值。
22max 1221max /1])1(5.0)1))(arctan(1[(a z p az v z a a z p z y x +-=+----==-σσ｛σ 式中，2max 32122212123/1/1/)1(/)1(83a F p d d E v E v F a π=+-+-=i i E v ,和)2,1(=j d j 分别为两个球的泊松比，弹性模量和直径。
用MA TLAB 表达形式写出上述公式并用下列给出的数值进行计算：)54.2(01.0)373(100,7.2,5.1,103,3.02172121厘米英寸，等于，克磅，等于in z lb F d d E E v v ====⨯====2、给出向量]19,7,51,37,29,0,5,47,3,17[----=x 。
编程使x 按照向量]5,17,19,29,51,47,37,7,3,0[----=y 的形式重新排列。
程序应使用于任意长度的向量。
对所有向量，数值0与负数放在一起。
也就是说，如果0为向量的一个元素，则它将作为y 向量的第一个元素。
3、如图所示曲柄滑块机构，滑块的位移为：22))sin(()cos(e a b a s --+=φφ位移s 是角度φ（以角度表示）的函数，当 3600,3.0,5.1,1≤≤===φe b a 时，使用plot(φ，s)函数绘制s 的图形。
4、由下式产生黄金分割数： ,......2,1,0],)251()251[(51=--+=n F n n n 使用fprintf 和disp 产生前16个数值，并在MA TLAB 命令窗口中显示如下内容：F0=0F1=1F2=1F3=2…….F15=6105、find 函数的实现给定由正数和负数构成的任意长度向量g,用程序实现下式：Indx=find(g>a)这里a 是由用户指定的一个值。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
力值。
｛σ x σ y
pmax [(1
z
a
arctan( ))(1
a
z
v1 )
0.5(1
z2 a2
)
1]
σz
p max 1 z2 / a 2
a
式中，
2
3F (1
3
v1 ) / E1
(1 v22 ) / E 2
第五题 a=4;
g=[4 4 7 10 -6 42 1 0]; % 输入函数 g 为一个行向量 indx=find(g>a)
第六题
t = 0:.01:2*pi; polar(t,2*ones(size(t)));
title('
极
坐
标
画
圆
')
z 0.01in(英寸，等于 2.54厘米 )
2 、 给 出 向 量 x [17, 3, 47,5,0,29, 37,51, 7,19] 。 编 程 使 x 按 照 向 量 y [ 0, 3, 7, 37, 47,51,29,19,17,5] 的形式重新排列。程序应使用于任意长度的向量。对所有向量，数
第四题 n=0:1:15; F= [n; (1/5^(1/2))*(((1+5^(1/2))/2).^n-((1-5^(1/2))/2).^n)]; fid = fopen('exp.txt','w'); fprintf(fid,'%2d F=%3.0f\n',F); fclose(fid); type exp.txt
值 0 与负数放在一起。也就是说，如果 0 为向量的一个元素，则它将作为 y 向量的第一个元素。
3、 如图所示曲柄滑块机构，滑块的位移为：
s a cos(φ) b2 (a sin(φ) e) 2
位移 s 是角度φ（以角度表示）的函数，当 plot( φ， s)函数绘制 s 的图形。
a 1,b 1.5, e 0.3,0 φ 360 时，使用
第二题 x=[17,-3,-47,5,0,29,-37,51,-7,19] indexzero=find(x==0);% x 是一个行向量 zero=x(indexzero); indexzheng=find(x>0); zheng=x(indexzheng); zhengp=-sort(-zheng);% 降序排列 indexfu=find(x<0); fu=x(indexfu); fup=-sort(-fu);% 降序排列 y=[zero,fup,zhengp]
4、 由下式产生黄金分割数：
Fn
11 [(
5)n
1 (
5 ) n ], n 0,1,2,......
52
2
使用 fprintf 和 disp 产生前 16 个数值，并在 MA TLAB 命令窗口中显示如下内容：
F0=0
F1=1
F2=1
F3=2
…… .
F15=610
5、 find 函数的实现 给定由正数和负数构成的任意长度向量
第三题 function s=zcj03(a,b,e) % 函数意义：曲柄滑块机构位移函数式 % 输入参数： a,,b,e % 输出参数： s a=0.5;b=1;e=0.2;% 该语句函数值可以在使用时，自己赋值 theta=-pi:0.1:pi; s=a.*cos(theta)+(b^2-(a.*sin(theta)-e).^2).^(1/2); plot(theta,s); end
8
1/ d1 1 / d2
p max
3F 2π a 2
vi , Ei 和 d j ( j 1,2) 分别为两个球的泊松比，弹性模量和直径。
用 MA TLAB 表达形式写出上述公式并用下列给出的数值进行计算：
v1 v2 0.3, E1 E2 3 107 , d1 1.5, d 2 2.7, F 100lb (磅，等于 373克)，
Indx=find(g>a) 这里 a 是由用户指定的一个值。本例中
g,用程序实现下式：
a=4，向量 g [ 4,4,7,10, 6,42,1,0] 。
6、 用极坐标画半径为 2 的圆 答案如下： 第一题 function [thetax,thetay,thetaz]= zcj01(v1,v2,E1,E2,d1,d2,F,z) % 函数功能：计算 xyz 方向应力 % 输入参数： v1,v2,E1,E2,d1,d2,F,z % 输出参数： thetax,thetay,thetaz v1=0.3;v2=0.3; E1=3*10^7;E2=3*10^7; d1=1.5;d2=2.7; F=37300; z=0.254; a=((3*F/8)*(((1-v1^2)/E1)+((1-v2^2)/E2))/(1/d1+1/d2))^(1/3); pmax=3*F/(2*pi*(a^2)); thetax=-pmax*((1-(z/a)*atan(a/z))*(1-v1)-0.5*(1+z^2/(a^2))^(-1)); thetay=-pmax*((1-(z/a)*atan(a/z))*(1-v1)-0.5*(1+z^2/(a^2))^(-1)); thetaz=-pmax/(1+z^2/(a^2)); end