matlab数学实验体会 |
您所在的位置:网站首页 › 数学实验课程心得体会 › matlab数学实验体会 |
matlab数学实验体会
时间:2024.3.11
Matlab学习心得 数学实验周我们学习了Matlab软件,这是一个十分实用和重要的软件。初次结识MATLAB,感觉这是一个很好的软件,语言简便,实用性强。作为一个新手,想要学习好这门语言,可以说还是比较难的。在我接触这门语言的这些天,一直在上面弄,除了会画几个简单的三维图形,其他的还是有待提高。在这个软件中,虽然有help。大家不要以为有了这个就万事大吉了,反而,从另一个方面也对我们大学生提出了两个要求――充实的课外基础和良好的英语基础。在现代,几乎所有好的软件都是来自国外,假如你不会外语,想学好是非常难的。 Matlab 语言是当今国际上科学界 (尤其是自动控制领域) 最具影响力、也是最有活力的软件。它起源于矩阵运算,并已经发展成一种高度集成的计算机语言。它提供了强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设计、便捷的与其他程序和语言接口的功能。 作为一种计算机语言,Matlab体现了与它价值的相符的优点: 1.编程简单使用方便。在这方面我感觉C语言也是一种简单的编程语言。只要入门就很好掌握,但是要学习一门语言不是那么容易的,到目前为止,可以说我还没入门,所以学习起这门语言来很吃力。相对C语言而言,Matlab的矩阵和向量操作功能是其他语言无法比拟的。在Matlab环境下,数组的操作与数的操作一样简单,基本数据单元是不需要指定维数的,不需要说明数据类型的矩阵,而其数学表达式和运算规则与通常的习惯相同。 2.函数库可任意扩充。由于Matlab语言库函数与用户文件的形式相同,用户文件可以像库函数一样随意调用,所以用户可任意扩充库函数。 3.语言简单内涵丰富。在此语言中,最重要的成分是函数,一般形式为:Function[a,b,c……]=fun(d,e,f……) Fun是自定义的函数名,只要不与库函数想重,并且符合字符串书写规则即可。 4、简便的绘图功能。MATLAB具有二维和三维绘图功能,使用方法简单。三维曲线是由plot3 (x,y,z)命令绘出的,看上去很简单的一个程序,相对C语言而言。极大的方便了绘图的工作和节省工作时间。 5. 丰富的工具箱。由于MATLAB 的开放性,许多领域的专家都为MATLAB 编写了各种程序工具箱。这些工具箱提供了用户在特别应用领域所需的许多函数,这使得用户不必花大量的时间编写程序就可以直接调用这些函数,达到事半功倍的效果。 在理论方面,在学习MATLAB过程中,我感觉到它和c语言有许多相似之处,他有c语言的特征,但是比c语言编程计算更加简单,适合于复杂的数学运算。但是MATLAB跟其他语言也有着很大的不同。 众所周知MATLAB是一个基于矩阵运算的软件,但是,真正在运用的时候,特别是在编程的时候,许多人往往没有注意到这个问题。在使用MATLAB时,受到了其他编程习惯的影响,特别是经常使用的C语言。因此,在MATLAB编程时,for循环(包括while循环)到处都是。这不仅是没有发挥MATLAB所长,还浪费了宝贵的时间。我这里想说的一点是,往往在初始化矩阵的时候注意到这个问题,懂得了使用矩阵而不是循环来赋值,但是,在其他环节上,就很容易疏忽,或者说,仍然没有摆脱C++、C的思想。多用help,see also, lookfor, get, set 等常用命令,尽量摆脱c编程的习惯,总爱用循环,能不用的循环的尽量不用,掌握矢量化的精髓。(1) help: 最有效的命令。其实,可以这样说吧,一遇 到什么问题,通常可以从 help 中找到答案。就先说说对help的一些常用方法。 1)命令窗口直接敲“help”,你就可以得到本地机器上matlab的基本的帮助信息。 2)对于某些不是很明确的命令,只知道大体所属范围,譬如说某个工具箱,直接在命令窗口中敲入 help toolboxname,一帮可以得到本工具箱有关的信息:版本号,函数名等。 (2) lookfor:可以说是 matlab中的google 。当我们很多什么头绪都没有的时候,我们可以求助于它,往往会收到意想不到的效果。譬如:曾经在gui编程的时候,遇到过这样一个问题:想拖动鼠标时,要出现一个方框,就像你在桌面上拖动鼠标,会出现虚线框一样。 当初我也刚开始一定都不知道该查找什么东西,后来想起用它了。于是, >> lookfor Rectangle。这样一条信息:GETRECT Select rectangle with mouse. get,set: GUI object 属性的帮手在GUI编程中,我们可能有时候想改变某些object的属性,或者想让它安装自己的想法实现,但是我们又不记得这些object的属性,更别提怎么设置他们的值了。这时,可以用 get(handles得到此对象的所有的属性及其当前值。用set(handles)可以得到对象所有可以设置的属性及其可能的取值找到我们需要的属性名字和可能的取值之后,就意义用 get(handles,‘propertyname’)取得此属性的值, 用set(handles,‘propertyname’,values)设置此对象此属性的值。 Edit: 查看m源文件的助手在应用matlab过程中,可能我们想看看它的m源文件,当然用editor定位打开也行,但是我经常采用的式直接在command窗口中用edit funname.m,就省去了定位的麻烦。 以上就是我学习MATLAB几个月以来的心得与体会,我自己感觉在理论方面自己理解的还是可以的,但是在实践中会经常遇到一些问题,而恰恰自己又束手无策。但是我经常上一些贴吧,那里有不少是使用MATLAB的高手,可以帮我解决不少问题,同时自己也学到了不少东西。 第二篇:matlab数学实验课后答案%Page20,ex1 (5) 等于[exp(1),exp(2);exp(3),exp(4)] (7) 3=1*3, 8=2*4 (8) a为各列最小值,b为最小值所在的行号 (10) 1>=4,false, 2>=3,false, 3>=2, ture, 4>=1,ture (11) 答案表明:编址第2元素满足不等式(30>=20)和编址第4元素满足不等式(40>=10) (12) 答案表明:编址第2行第1列元素满足不等式(30>=20)和编址第2行第2列元素满足不等式(40>=10) %Page20, ex2 (1)a, b, c的值尽管都是1, 但数据类型分别为数值,字符, 逻辑, 注意a与c相等, 但他们不等于b (2)double(fun)输出的分别是字符a,b,s,(,x,)的ASCII码 %Page20,ex3 >> r=2;p=0.5;n=12; >> T=log(r)/n/log(1+0.01*p) T = 11.5813 %Page20,ex4 >> x=-2:0.05:2;f=x.^4-2.^x; >> [fmin,min_index]=min(f) fmin = -1.3907 %最小值 min_index = 54 %最小值点编址 >> x(min_index) ans = 0.6500 %最小值点 >> [f1,x1_index]=min(abs(f)) %求近似根--绝对值最小的点 f1 = 0.0328 x1_index = 24 >> x(x1_index) ans = -0.8500 >> x(x1_index)=[];f=x.^4-2.^x; %删去绝对值最小的点以求函数绝对值次小的点 >> [f2,x2_index]=min(abs(f)) %求另一近似根--函数绝对值次小的点 f2 = 0.0630 x2_index = 65 >> x(x2_index) ans = 1.2500 %Page20,ex5 >> z=magic(10) z = 92 99 1 8 15 67 74 51 58 40 98 80 7 14 16 73 55 57 64 41 4 81 88 20 22 54 56 63 70 47 85 87 19 21 3 60 62 69 71 28 86 93 25 2 9 61 68 75 52 34 17 24 76 83 90 42 49 26 33 65 23 5 82 89 91 48 30 32 39 66 79 6 13 95 97 29 31 38 45 72 10 12 94 96 78 35 37 44 46 53 11 18 100 77 84 36 43 50 27 59 >> sum(z) ans = 505 505 505 505 505 505 505 505 505 505 >> sum(diag(z)) ans = 505 >> z(:,2)/sqrt(3) ans = 57.1577 46.1880 46.7654 50.2295 53.6936 13.8564 2.8868 3.4641 6.9282 10.3923 >> z(8,:)=z(8,:)+z(3,:) z = 92 99 1 8 15 67 74 51 58 40 98 80 7 14 16 73 55 57 64 41 4 81 88 20 22 54 56 63 70 47 85 87 19 21 3 60 62 69 71 28 86 93 25 2 9 61 68 75 52 34 17 24 76 83 90 42 49 26 33 65 23 5 82 89 91 48 30 32 39 66 83 87 101 115 119 83 87 101 115 119 10 12 94 96 78 35 37 44 46 53 11 18 100 77 84 36 43 50 27 59 %Page 40 ex1 先在编辑器窗口写下列M函数,保存为eg2_1.m function [xbar,s]=ex2_1(x) n=length(x); xbar=sum(x)/n; s=sqrt((sum(x.^2)-n*xbar^2)/(n-1)); 例如 >>x=[81 70 65 51 76 66 90 87 61 77]; >>[xbar,s]=ex2_1(x) xbar = 72.4000 s = 12.1124 %Page 40 ex2 s=log(1);n=0; while se k=k+1;F(k)=F(k-1)+F(k-2); x=F(k)/F(k-1); end a,x,k 计算至k=21可满足精度 %Page 40 ex4 clear;tic;s=0; for i=1:1000000 s=s+sqrt(3)/2^i; end s,toc tic;s=0;i=1; while i1.1)+x.*(x=-1.1)-1.1*(x1); p=p+b*exp(-y.^2-6*x.^2).*(x+y>-1).*(x+y A=[1 2 3;4 5 6;7 8 0];C=[2 -5 -22;-5 -24 -56;-22 -56 -16]; >> X=lyap(A,C) X = 1.0000 -1.0000 -0.0000 -1.0000 2.0000 1.0000 -0.0000 1.0000 7.0000 %Chapter 3 %Exercise 1 >> a=[1,2,3];b=[2,4,3];a./b,a.\b,a/b,a\b ans = 0.5000 0.5000 1.0000 ans = 2 2 1 ans = 0.6552 %一元方程组x[2,4,3]=[1,2,3]的近似解 ans = 0 0 0 0 0 0 0.6667 1.3333 1.0000 %矩阵方程[1,2,3][x11,x12,x13;x21,x22,x23;x31,x32,x33]=[2,4,3]的特解 %Exercise 2(1) >> A=[4 1 -1;3 2 -6;1 -5 3];b=[9;-2;1]; >> rank(A), rank([A,b]) %[A,b]为增广矩阵 ans = 3 ans = 3 %可见方程组唯一解 >> x=A\b x = 2.3830 1.4894 2.0213 %Exercise 2(2) >> A=[4 -3 3;3 2 -6;1 -5 3];b=[-1;-2;1]; >> rank(A), rank([A,b]) ans = 3 ans = 3 %可见方程组唯一解 >> x=A\b x = -0.4706 -0.2941 %Exercise 2(3) >> A=[4 1;3 2;1 -5];b=[1;1;1]; >> rank(A), rank([A,b]) ans = 2 ans = 3 %可见方程组无解 >> x=A\b x = 0.3311 -0.1219 %最小二乘近似解 %Exercise 2(4) >> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1 2 3]';%注意b的写法 >> rank(a),rank([a,b]) ans = 3 ans = 3 %rank(a)==rank([a,b])> a\b ans = 1 1 0 %一个特解 %Exercise 3 >> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1,2,3]'; >> x=null(a),x0=a\b x = -0.6255 0.6255 -0.2085 0.4170 x0 = 1 1 %通解kx+x0 另一种解法 >> rref([a,b]) ans = 1.0000 0 0 1.5000 1.0000 0 1.0000 0 -1.5000 0 0 0 1.0000 0.5000 1.0000 %通解k[1.5,-1.5,0.5,1]'+[1,0,1,0]' %Exercise 4 %一年后城市=上年城市×99%+上年农村×5%;一年后农村=上年城市×1%+上年农村×95% >> x0=[0.2 0.8]';a=[0.99 0.05;0.01 0.95]; >> x1=a*x, x2=a^2*x, x10=a^10*x >> x=x0;for i=1:1000,x=a*x;end,x x = 0.8333 0.1667 >> x0=[0.7 0.3]'; >> x=x0;for i=1:1000,x=a*x;end,x x = 0.8333 0.1667 %注意结果不受x0改变的影响 >> [v,e]=eig(a) v = 0.9806 -0.7071 0.1961 0.7071 e = 1.0000 0 0 0.9400 >> v(:,1)./x ans = 1.1767 1.1767 %v(:,1)与x成比例,说明x是最大特征值对应的特征向量 %Exercise 5 %用到公式(3.11)(3.12) >> B=[6,2,1;2.25,1,0.2;3,0.2,1.8];x=[25 5 20]'; >> C=B/diag(x) C = 0.2400 0.4000 0.0500 0.0900 0.2000 0.0100 0.1200 0.0400 0.0900 >> A=eye(3,3)-C A = 0.7600 -0.4000 -0.0500 -0.0900 0.8000 -0.0100 -0.1200 -0.0400 0.9100 >> D=[17 17 17]';x=A\D x = 37.5696 25.7862 24.7690 %Exercise 6(1) >> a=[4 1 -1;3 2 -6;1 -5 3];det(a),inv(a),[v,d]=eig(a) ans = -94 ans = 0.2553 -0.0213 0.0426 0.1596 -0.1383 -0.2234 0.1809 -0.2234 -0.0532 v = 0.0185 -0.9009 -0.3066 -0.7693 -0.1240 -0.7248 -0.6386 -0.4158 0.6170 d = -3.0527 0 0 0 3.6760 0 0 0 8.3766 %Exercise 6(2) >> a=[1 1 -1;0 2 -1;-1 2 0];det(a),inv(a),[v,d]=eig(a) ans = 1 ans = 2.0000 -2.0000 1.0000 1.0000 -1.0000 1.0000 2.0000 -3.0000 2.0000 v = -0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i -0.5773 0.5774 0.5774 -0.5774 0.5773 - 0.0000i 0.5773 + 0.0000i d = 1.0000 0 0 0 1.0000 + 0.0000i 0 0 0 1.0000 - 0.0000i %Exercise 6(3) >> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10] A = 5 7 6 5 7 10 8 7 6 8 10 9 5 7 9 10 >> det(A),inv(A), [v,d]=eig(A) ans = 1 ans = 68.0000 -41.0000 -17.0000 10.0000 -41.0000 25.0000 10.0000 -6.0000 -17.0000 10.0000 5.0000 -3.0000 10.0000 -6.0000 -3.0000 2.0000 v = 0.8304 0.0933 0.3963 0.3803 -0.5016 -0.3017 0.6149 0.5286 -0.2086 0.7603 -0.2716 0.5520 0.1237 -0.5676 -0.6254 0.5209 d = 0.0102 0 0 0 0 0.8431 0 0 0 0 3.8581 0 0 0 0 30.2887 %Exercise 6(4)(以n=5为例,n=其他值类似) %关键是矩阵的定义,当n很大时不能直接写出。 %方法一(三个for) n=5; for i=1:n, a(i,i)=5;end for i=1:(n-1),a(i,i+1)=6;end for i=1:(n-1),a(i+1,i)=1;end a %方法二(一个for) n=5;a=zeros(n,n); a(1,1:2)=[5 6]; for i=2:(n-1),a(i,[i-1,i,i+1])=[1 5 6];end a(n,[n-1 n])=[1 5]; a %方法三(不用for) n=5;a=diag(5*ones(n,1)); b=diag(6*ones(n-1,1)); c=diag(ones(n-1,1)); a=a+[zeros(n-1,1),b;zeros(1,n)]+[zeros(1,n);c,zeros(n-1,1)] %下列计算 >> det(a) ans = 665 >> inv(a) ans = 0.3173 -0.5865 1.0286 -1.6241 1.9489 -0.0977 0.4887 -0.8571 1.3534 -1.6241 0.0286 -0.1429 0.5429 -0.8571 1.0286 -0.0075 0.0376 -0.1429 0.4887 -0.5865 0.0015 -0.0075 0.0286 -0.0977 0.3173 >> [v,d]=eig(a) v = -0.7843 -0.7843 -0.9237 0.9860 -0.9237 0.5546 -0.5546 -0.3771 -0.0000 0.3771 -0.2614 -0.2614 0.0000 -0.1643 0.0000 0.0924 -0.0924 0.0628 -0.0000 -0.0628 -0.0218 -0.0218 0.0257 0.0274 0.0257 d = 0.7574 0 0 0 0 0 9.2426 0 0 0 0 0 7.4495 0 0 0 0 0 5.0000 0 0 0 0 0 2.5505 %Exercise 7(1) >> a=[4 1 -1;3 2 -6;1 -5 3];[v,d]=eig(a) v = 0.0185 -0.9009 -0.3066 -0.7693 -0.1240 -0.7248 -0.6386 -0.4158 0.6170 d = -3.0527 0 0 0 3.6760 0 0 0 8.3766 >> det(v) ans = -0.9255 %v行列式正常, 特征向量线性相关,可对角化 >> inv(v)*a*v %验算 ans = -3.0527 0.0000 -0.0000 0.0000 3.6760 -0.0000 -0.0000 -0.0000 8.3766 >> [v2,d2]=jordan(a) %另一方法:用jordan v2 = 0.0798 0.0076 0.9127 0.1886 -0.3141 0.1256 -0.1605 -0.2607 0.4213 %用jordan特征向量不同 d2 = 8.3766 0 0 0 -3.0527 - 0.0000i 0 0 0 3.6760 + 0.0000i >> v2\a*v2 ans = 8.3766 0 0.0000 0.0000 -3.0527 0.0000 0.0000 0.0000 3.6760 >> v(:,1)./v2(:,2) %对应相同特征值的特征向量成比例,说明两种方法结果本质上相同 ans = 2.4491 2.4491 2.4491 %Exercise 7(2) >> a=[1 1 -1;0 2 -1;-1 2 0];[v,d]=eig(a) v = -0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i -0.5773 0.5774 0.5774 -0.5774 0.5773 - 0.0000i 0.5773 + 0.0000i d = 1.0000 0 0 0 1.0000 + 0.0000i 0 0 0 1.0000 - 0.0000i >> det(v) ans = -5.0566e-028 -5.1918e-017i %v的行列式接近0, 特征向量线性相关,不可对角化 >> [v,d]=jordan(a) v = 1 0 1 1 0 0 1 -1 0 d = 1 1 0 0 1 1 0 0 1 %jordan标准形不是对角的,所以不可对角化 %Exercise 7(3) >> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10] A = 5 7 6 5 7 10 8 7 6 8 10 9 5 7 9 10 >> [v,d]=eig(A) v = 0.8304 0.0933 0.3963 0.3803 -0.5016 -0.3017 0.6149 0.5286 -0.2086 0.7603 -0.2716 0.5520 0.1237 -0.5676 -0.6254 0.5209 d = 0.0102 0 0 0 0 0.8431 0 0 0 0 3.8581 0 0 0 0 30.2887 >> inv(v)*A*v ans = 0.0102 0.0000 -0.0000 0.0000 0.0000 0.8431 -0.0000 -0.0000 -0.0000 0.0000 3.8581 -0.0000 -0.0000 -0.0000 0 30.2887 %但本题用jordan不行, 原因未知 %Exercise 7(4)参考6(4)和7(1), 略 %Exercise 8 只有(3)对称且特征值全部大于零, 所以只有(3)是正定矩阵. %Exercise 9 >> a=[4 -3 1 3;2 -1 3 5;1 -1 -1 -1;3 -2 3 4;7 -6 -7 0] >> rank(a) ans = 3 >> rank(a(1:3,:)) ans = 2 >> rank(a([1 2 4],:)) %1,2,4行为最大无关组 ans = 3 >> b=a([1 2 4],:)';c=a([3 5],:)'; >> b\c %线性表示的系数 ans = 0.5000 5.0000 -0.5000 1.0000 0 -5.0000 %Exercise 10 >> a=[1 -2 2;-2 -2 4;2 4 -2] >> [v,d]=eig(a) v = 0.3333 0.9339 -0.1293 0.6667 -0.3304 -0.6681 -0.6667 0.1365 -0.7327 d = -7.0000 0 0 0 2.0000 0 0 0 2.0000 >> v'*v ans = 1.0000 0.0000 0.0000 0.0000 1.0000 0 0.0000 0 1.0000 %可见v确实是正交矩阵 %Exercise 11 %设经过6个电阻的电流分别为i1, ..., i6. 列方程组如下 %20-2i1=a; 5-3i2=c; a-3i3=c; a-4i4=b; c-5i5=b; b-3i6=0; %i1=i3+i4;i5=i2+i3;i6=i4+i5; %计算如下 >> A=[1 0 0 2 0 0 0 0 0;0 0 1 0 3 0 0 0 0;1 0 -1 0 0 -3 0 0 0; 1 -1 0 0 0 0 -4 0 0; 0 -1 1 0 0 0 0 -5 0;0 1 0 0 0 0 0 0 -3; 0 0 0 1 0 -1 -1 0 0;0 0 0 0 -1 -1 0 1 0; 0 0 0 0 0 0 -1 -1 1]; >>b=[20 5 0 0 0 0 0 0 0]'; A\b ans = 13.3453 6.4401 8.5420 3.3274 -1.1807 1.6011 1.7263 0.4204 2.1467 %Exercise 12 >> A=[1 2 3;4 5 6;7 8 0]; >> left=sum(eig(A)), right=sum(trace(A)) left = 6.0000 right = 6 >> left=prod(eig(A)), right=det(A) %原题有错, (-1)^n应删去 left = 27.0000 right = 27 >> fA=(A-p(1)*eye(3,3))*(A-p(2)*eye(3,3))*(A-p(3)*eye(3,3)) fA = 1.0e-012 * 0.0853 0.1421 0.0284 0.1421 0.1421 0 -0.0568 -0.1137 0.1705 >> norm(fA) ans = 2.9536e-013 %f(A)范数接近0,也就是说fA是0矩阵。 %Exercise 1(1) roots([1 1 1]) %Exercise 1(2) roots([3 0 -4 0 2 -1]) %Exercise 1(3) p=zeros(1,24); p([1 17 18 22])=[5 -6 8 -5]; %这样比直接写24个系数简短。 roots(p) %Exercise 1(4) p1=[2 3]; p2=conv(p1, p1); p3=conv(p1, p2); p3(end)=p3(end)-4; %原p3最后一个分量-4 roots(p3) %Exercise 2 fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x'); fzero(fun,2) %Exercise 3 fun=inline('x^4-2^x'); fplot(fun,[-2 2]);grid on; fzero(fun,-1),fzero(fun,1),fminbnd(fun,0.5,1.5) %Exercise 4 fun=inline('x*sin(1/x)','x'); fplot(fun, [-0.1 0.1]); x=zeros(1,10);for i=1:10, x(i)=fzero(fun,(i-0.5)*0.01);end; x=[x,-x] %Exercise 5 fun=inline('[9*x(1)^2+36*x(2)^2+4*x(3)^2-36;x(1)^2-2*x(2)^2-20*x(3);16*x(1)-x(1)^3-2*x(2)^2-16*x(3)^2]','x'); [a,b,c]=fsolve(fun,[0 0 0]) %Exercise 6 fun=@(x)[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))]; [a,b,c]=fsolve(fun,[0.5 0.5]) %Exercise 7 clear; close; t=0:pi/100:2*pi; x1=2+sqrt(5)*cos(t); y1=3-2*x1+sqrt(5)*sin(t); x2=3+sqrt(2)*cos(t); y2=6*sin(t); plot(x1,y1,x2,y2); grid on; %作图发现4个解的大致位置,然后分别求解 y1=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.5,2]) y2=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.8,-2]) y3=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[3.5,-5]) y4=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[4,-4]) %Exercise 8(1) clear; str='x.^2.*sin(x.^2-x-2)'; %注意数组点运算 fun=inline(str); fplot(fun,[-2 2]);grid on; %作图观察 x(1)=-2; x(3)=fminbnd(fun,-1,-0.5); x(5)=fminbnd(fun,1,2); fun2=inline(['-',str]); x(2)=fminbnd(fun2,-2,-1); x(4)=fminbnd(fun2,-0.5,0.5); x(6)=2 feval(fun,x) %答案: 以上x(1)(3)(5)是局部极小,x(2)(4)(6)是局部极大,从最后一句知道x(1)全局最小, x(2)最大。 %Exercise 8(2) clear; fun=inline('3*x.^5-20*x.^3+10'); fplot(fun,[-3 3]);grid on;%作图观察 x(1)=-3; x(3)=fminsearch(fun,2.5); fun2=inline('-(3*x.^5-20*x.^3+10)'); x(2)=fminsearch(fun2,-2.5); x(4)=3; feval(fun,x) %Exercise 8(3) fun=inline('abs(x^3-x^2-x-2)'); fplot(fun,[0 3]);grid on;%作图观察 fminbnd(fun,1.5,2.5) fun2=inline('-abs(x^3-x^2-x-2)'); fminbnd(fun2,0.5,1.5) %Exercise 9 close; x=-2:0.1:1;y=-7:0.1:1; [x,y]=meshgrid(x,y); z=y.^3/9+3*x.^2.*y+9*x.^2+y.^2+x.*y+9; mesh(x,y,z);grid on;%作图观察, 可看到[0 0]附近极小值,[0 -5]附近极大值 fun=inline('x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9'); x=fminsearch(fun,[0 0])%求极小值 fun2=inline('-(x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9)'); x=fminsearch(fun2,[0 -5])%求极大值 %Exercise 10 clear;t=0:24; c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ... 31 32 31 29 27 25 24 22 20 18 17 16]; %第(1)小题 p2=polyfit(t,c,2) %第(2)小题 p3=polyfit(t,c,3) %第(3)小题 fun=inline('a(1)*exp(a(2)*(t-14).^2)','a','t'); a=lsqcurvefit(fun,[0 0],t,c)%初值可以试探,这里用a=0, b=0试探 f=feval(fun, a,t) norm(f-c)%拟合效果 plot(t,c,t,f) %作图检验 %第(4)小题 fun2=inline('b(1)*sin(pi/12*t+b(2))+20','b','t');%原题修改为f(x)+20,结果较合理 b=lsqcurvefit(fun2,[0 0],t,c) figure f2=feval(fun2, b,t) norm(f2-c)%拟合效果 plot(t,c,t,f2) %作图检验 %Exercise 11 fun=inline('(1-x)*sqrt(10.52+x)-3.06*x*sqrt(1+x)*sqrt(5)'); x=fzero(fun, 0, 1) %Exercise 12 r=5.04/12/100;N=20*12; x=7500*180 %房屋总价格 y=x*0.3 %首付款额 x0=x-y%贷款总额 a=(1+r)^N*r*x0/((1+r)^N-1)%月付还款额 %以下组合贷款 r1=4.05/12/100;x1=10*10000;%公积金贷款 a1=(1+r1)^N*r1*x1/((1+r1)^N-1) x2=x0-x1%商业贷款 a2=(1+r)^N*r*x2/((1+r)^N-1) a=a1+a2 %Exercise 13 %列方程th*R^2+(pi-2*th)*r^2-R*r*sin(th)=pi*r^2/2 %化简得sin(2*th)-2*th*cos(2*th)=pi/2 %以下Matlab计算 clear;fun= inline('sin(2*th)-2*th*cos(2*th)-pi/2','th') th=fsolve(fun,pi/4) R=20*cos(th) %Exercise 14 %先在Editor窗口写M函数保存 function x=secant(fname,x0,x1,e) while abs(x0-x1)>e, x=x1-(x1-x0)*feval(fname,x1)/(feval(fname,x1)-feval(fname,x0)); x0=x1;x1=x; end %再在指令窗口 fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x'); secant(fun,1,2,1e-8) %Exercise 15 %作系数为a,初值为xo,从第m步到第n步迭代过程的M函数: function f=ex4_15fun(a,x0,m,n) x(1)=x0; y(1)=a*x(1)+1;x(2)=y(1); if mm, plot([x(i),x(i),x(i+1)],[y(i-1),y(i),y(i)]); end end hold off; %M脚本文件 subplot(2,2,1);ex4_15fun(0.9,1,1,20); subplot(2,2,2);ex4_15fun(-0.9,1,1,20); subplot(2,2,3);ex4_15fun(1.1,1,1,20); subplot(2,2,4);ex4_15fun(-1.1,1,1,20); %Exercise 16 %设夹角t, 问题转化为 min f=5/sin(t)+10/cos(t) %取初始值pi/4, 计算如下 fun=@(t)5/sin(t)+10/cos(t); [t,f]=fminsearch(fun, pi/4) t = 0.6709 f = 20.8097 %Exercise 17 %提示:x(k+2)=f(x(k))=a^2*x(k)*(1-x(k))*(1-a*x(k)*(1-x(k))) %计算平衡点x %|f'(x)|> ex4_18(0.9,1,20) >> ex4_18(-0.9,1,20) >> ex4_18(1.1,1,20) >> ex4_18(-1.1,1,20) %Exercise 19 clear; close; x(1)=0; y(1)=0; for k=1:3000 x(k+1)=1+y(k)-1.4*x(k)^2; y(k+1)=0.3*x(k); end plot(x(1000:1500),y(1000:1500),'+g');hold on plot(x(1501:2000),y(1501:2000),'.b'); plot(x(2001:2500),y(2001:2500),'*y'); plot(x(2501:3001),y(2501:3001),'.r'); %Exercise 1 x=[0 4 10 12 15 22 28 34 40]; y=[0 1 3 6 8 9 5 3 0]; trapz(x,y) %Exercise 2 x=[0 4 10 12 15 22 28 34 40]; y=[0 1 3 6 8 9 5 3 0]; diff(y)./diff(x) %Exercise 3 xa=-1:0.1:1;ya=0:0.1:2; [x,y]=meshgrid(xa,ya); z=x.*exp(-x.^2 -y.^3); [px,py] = gradient(z,xa,ya); px %Exercise 4 t=0:0.01:1.5; x=log(cos(t)); y=cos(t)-t.*sin(t); dydx=gradient(y,x) plot(x,dydx) %dydx函数图,作图观察x=-1时,dydx的值约0.9 %以下是更精确的编程计算方法 [x_1,id]=min(abs(x-(-1)));%找最接近x=-1的点,id为这个点的下标 dydx(id) %Exercise 5(1) fun=@(x)1/sqrt(2*pi)*exp(-x.^2/2); quadl(fun,0,1) 或用trapz x=linspace(0,1,100); y=1/sqrt(2*pi)*exp(-x.^2/2); trapz(x,y) %Exercise 5(2) fun=inline('exp(2*x).*cos(x).^3'); quadl(fun,0,2*pi) 或用trapz x=linspace(0,2*pi,100); y=exp(2*x).*cos(x).^3; trapz(x,y) %Exercise 5(3) fun=@(x)x.*log(x.^4).*asin(1./x.^2); quadl(fun,1,3) 或用trapz x=1:0.01:3; y=feval(fun,x); trapz(x,y) %Exercise 5(4) fun=@(x)sin(x)./x; quadl(fun,1e-10,1) %注意由于下限为0,被积函数没有意义,用很小的1e-10代替 %Exercise 5(5) fun=inline('x.^(-x)','x'); quadl(fun,1e-10,1) %注意由于下限为0,被积函数没有意义,用很小的1e-10代替 %Exercise 5(6) fun=inline('sqrt(1+r.^2.*sin(th))','r','th'); dblquad(fun,0,1,0,2*pi) %Exercise 5(7) %先在Editer窗口建立90页函数dblquad2,再在Command窗口 clear; fun=@(x,y)1+x+y.^2; clo=@(x)-sqrt(2*x-x.^2); dhi=@(x)sqrt(2*x-x.^2); dblquad2(fun,0,2,clo,dhi,100) %Exercise 6 t=linspace(0,2*pi,100); x=2*cos(t);y=3*sin(t); dx=gradient(x,t); dy=gradient(y,t); f=sqrt(dx.^2+dy.^2); trapz(t,f) %Exercise 7 %先算出z的梯度 dz/dx=(1-2*x^2)*exp(-x^2-y^2),dz/dy=(1-2*y^2)*exp(-x^2-y^2); %根据曲面面积公式 fun=@(x,y)sqrt(1+((1-2*x.^2).*exp(-x.^2-y.^2)).^2+((1-2*y.^2)*exp(-x.^2-y.^2)).^2); dblquad(fun,-1,1,0,2) %或者用下列纯粹离散化解法 xa=linspace(-1,1);ya=linspace(0,2); [x,y]=meshgrid(xa,ya); z=x.*exp(-x.^2-y.^2); [zx,zy]=gradient(z,xa,ya); f=sqrt(1+zx.^2+zy.^2); s=0; for i=2:length(xa) for j=2:length(ya) s=s+(xa(i)-xa(i-1))*(ya(j)-ya(j-1))*(f(i,j)+f(i-1,j)+f(i,j-1)+f(i-1,j-1))/4;%每个近似长方体高用四顶点平均值 end end s %Exercise 8 funl=inline('-(-x).^0.2.*cos(x)'); funr=inline('x.^0.2.*cos(x)'); quadl(funl,-1,0)+quadl(funr,0,1) %Exercise 9 (以I32为例) fun=@(x)abs(sin(x)); h=0.1;x=0:h:32*pi;y=feval(fun,x);t1=trapz(x,y) h=pi;x=0:h:32*pi;y=feval(fun,x);t2=trapz(x,y)%步长与周期一致,结果失真 q1=quad(fun,0,32*pi) q2=quadl(fun,0,32*pi) %Exercise 10(1) 先在Editer窗口建立88页函数deriv,再在Command窗口 fun=inline('x.^2.*sin(x.^2+3*x-4)','x'); deriv(fun,[1.3 1.5],0.1,1e-3) %取0.1为初始步长 %注:书后习题答案错,1.3处导数应为2.4177,1.5处导数应为-11.3330 %Exercise 10(2) %先在程序编辑器,写下列函数,保存为ex5_10_2f function d=ex5_10_2f(fname,a,h0,e) h=h0;d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h); d0=d+2*e; while abs(d-d0)>e d0=d;h0=h;h=h0/2; d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h); end %再在指令窗口执行 fun=inline('x.^2*sin(x.^2-x-2)','x'); d=ex5_10_2f(fun,1.4,0.1,1e-3) %Exercise 11 %提示:f上升时,f'>0;f下降时,f'0) %解得户外温度c=33,比例系数k=0.05. %Excercise 8 %设市场占有率x, 增速x', 相对增长率x'/x. 建立微分方程 x'/x=0.5*(1-x),x(0)=0.05 %解析解x(t)=1/(1+19*exp(-1/2*t)),可手工求或用Ch7方法 dsolve('Dx=0.5*x*(1-x)','x(0)=0.05','t') %以下数值解 clear; fun=inline('0.5*(1-x)*x','t','x'); [t,x]=ode45(fun,[0 10],0.05); %作图比较 close; fplot('1/(1+19*exp(-1/2*t))',[0,10]);hold on; %实线为解析解 plot(t,x,'o');hold off; %圆圈为数值解,发现与解析解基本重合 %找市场占有率达到80%的时间可从图上观察,也可用以下方法 id=min(find(x>0.8)); t(id) %Excercise 9 %微分方程组 V'(t)=K(t)*V(t)^a,K'(t)=-b*K(t) %令y(1)=V,y(2)=K. 编程计算如下 clear;close; a=1;b=0.1; dy=@(t,y)[y(2)*y(1)^a;-b*y(2)]; [t,y]=ode45(dy,[0 150],[1;2]) plot(t,y(:,1)) %从图中可以看出有一条渐近线,约为4.8567*10^8. %时间区间右端点tf=150是试出来的,通过tf=10,20,50,100,150逐步增大看出趋势为止 %第(3)小题可从图上观测t大约为30 figure;dv=y(:,1).*y(:,2).^a;plot(t,dv);grid on; %也可如下计算 [mv,id]=max(dv);t(id) %计算结果30 %第(2)小题理论上也可从图上观测,但是由于t范围太大,精度不高。 %为了提高精度,缩小t范围重新计算。 [t,y]=ode45(dy,[0 1],[1;2]) id=min(find(y(:,1)>2)); t(id) %计算结果0.375 %第(4)小题,只要将上述a改为2/3计算,相应t右端点tf可取80. 计算结果3个小题答案依次451, 0.4, 15. %本题求解也可用Ch7解析方法 [K,V]=dsolve('DK=-0.1*K','DV=K*V','K(0)=2','V(0)=1','t') %得到V=exp(2)^10*exp(-20*exp(-1/10*t)),据此可以得到更精确的结果 syms t;limit(V,t,inf) Vfun=inline([char(V),'-2'],'t'); fzero(Vfun,1) V2=diff(V,t,2); Vdfun=inline(char(V2),'t'); fzero(Vdfun,20) %Excercise 10 %先在Editor窗口写M-函数 function dx=lorenfun(t,x) sig=10;r=28;b=8/3; dx=[-sig sig 0;r -1 -x(1);0 x(1) -b]*x; %再在Commond窗口执行 clear;close; [t,x]=ode45(@lorenfun,[0,20],[10;-10;-10]); [t(end),x(end,:)] plot3(x(:,1),x(:,2),x(:,3),'b');hold on;%蓝色线 [t,x]=ode45(@lorenfun,[0,20],[10.001;-10;-10]); [t(end),x(end,:)] plot3(x(:,1),x(:,2),x(:,3),'g');%绿色 [t,x]=ode45(@lorenfun,[0,20],[-10;-10;-10]); [t(end),x(end,:)] plot3(x(:,1),x(:,2),x(:,3),'r');hold off;%红色 %Excercise 11 %设电量为Q(t), 电流为 I(t), 根据基尔霍夫第二定律, % L*I'(t)+R*I(t)+Q(t)/C=E % 因为I(t)=Q'(t), 对上式求导得 % I''(t)+R/L*I'(t)+I(t)/L/C=0 %那么电压u(t)=R*I(t)满足 % u''(t)+R/L*u'(t)+u(t)/L/C=0 % 初始u(0)=20, u'(0)=0; % 现在用Matlab求解,令y(1)=u, y(2)=u', 那么 y'(1)=y(2),y'(2)=-R/L*y(2)-y(1)/L/C; clear;close; R=100;L=0.1;C=0.2; fun=@(t,y)[y(2);-R/L*y(2)-y(1)/L/C]; [t,u]=ode45(fun,[0 50],[20;0]); %时间终点50是尝试性的, 势即可 [t,u(:,1)] %第一列是电压 plot(t,u(:,1)) 看出趋于0的趋 %Excercise 12 %先在Editor窗口写M-函数 function dx=volterafun(t,x) a1=3;b1=2;a2=2.5;b2=1; dx(1)=x(1)*(a1-b1*x(2)); dx(2)=-x(2)*(a2-b2*x(1)); dx=dx(:); %再在Commond窗口执行 clear;close; [t,x]=ode45(@volterafun,[0,10],[1;1]); %时间终点10是尝试性的,看出2个以上周期即可 subplot(1,2,1);plot(t,x(:,1),t,x(:,2)); %解曲线 subplot(1,2,2);plot(x(:,1),x(:,2)); %相轨线 %Excercise 13 %先在Editor窗口写M-函数 function [x,y]=nark4(dyfun,xspan,y0,h) %用途: 4阶经典Runge-Kutta法解常微分方程 y'=f(x,y),y(x0)=y0 %格式: [x,y]=nark4(dyfun,xspan,y0,h),dyfun为函数f(x,y),xspan为求解区间[x0,xN] % y0为初值y(x0),h为步长,x返回节点,y返回数值解 x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2); k4=feval(dyfun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6; end x=x';y=y'; %再在Commond窗口执行 fun=inline('x+y','x','y'); [t,y]=nark4(fun,[0 3],1,0.1); %步长h取为0.1 [t,y] %Excercise 14 %令y(1)=x,y(2)=x', 化为一阶方程组 clear;close; g=9.8;L=150;m=70;k=10;c0=0.00324; fun=@(t,y)[y(2);-c0*sign(y(2))*y(2)^2-k/m*(y(1)-L)*(y(1)>L)+g]; [t,x]=ode45(fun,[0 5.988 11.18 11.47],[0;0]);[t,x] d2x=-c0*sign(x(3,2))*x(3,2)^2-k/m*(x(3,1)-L)*(x(3,1)>L)+g [t,x]=ode45(fun,[0 30],[0;0]); d2x=-c0*sign(x(:,2)).*x(:,2).^2-k/m*(x(:,1)-L).*(x(:,1)>L)+g; subplot(1,3,1);plot(t,x(:,1));title('位移') subplot(1,3,2);plot(t,x(:,2));title('速度') subplot(1,3,3);plot(t,d2x);title('加速度') max(x(:,1)), min(d2x) %Chapter 7 %Exercise 1 syms ph th; a=sin(ph)*cos(th)-cos(ph)*sin(th)-sin(ph-th); simple(a) %化简后差的结果为0 %Exercise 2 syms x;s=x^4-5*x^3+5*x^2+5*x-6; factor(s) %Exercise 3 syms a;A=[1 2;2 a]; iA=inv(A),[v,d]=eig(A) %Exercise 4 syms x y; limit((3^x+9^x)^(1/x),x,inf) s1=limit(x*y/(sqrt(x*y+1)-1),x,0);s2=limit(s1,y,0) %Exercise 5 syms k n x;s1=symsum(k^2,k,1,n);s1=simple(s1) s2=symsum(k^(-2),k,1,inf);s2=simple(s2) s3=symsum(1/(2*n+1)/(2*x+1)^(2*n+1),n,0,inf);s3=simple(s3) %Exercise 6 syms x y z;s=sin(x^2*y*z); s=diff(s,x,2); s=diff(s,y,1); s=subs(s,{x,y,z},{1,1,3}) %Exercise 7 syms x; s=exp(x);taylor(s,8,0,x) s=log(1+x);taylor(s,8,0,x) s=sin(x);taylor(s,8,0,x) s=log(x+sqrt(1+x^2));taylor(s,8,0,x) %Exercise 8 %第四章习题8(1) %先用符号运算求偏导数 syms x;f=x^2*sin(x^2-x-2); fx=diff(f,x) %根据计算结果得方程组. 求解方程组 sx=solve(fx,x) %得 0 和 -0.7315。考虑2阶导数 fx2=diff(fx,x),vpa(subs(fx2,sx)) %知道0是极大值,-0.7315是极小值 %注意:解析方法并没有求得所有极值点 %第四章习题8(2)类似 %第四章习题8(3)因为有绝对值(不可导),只能去掉绝对值来用diff分析。 %以下第四章习题9 %先用符号运算求偏导数 syms x y;f=y^3/9+3*x^2*y+9*x^2+y^2+x*y+9; fx=diff(f,x),fy=diff(f,y) %根据计算结果得方程组. 求解方程组 [sx,sy]=solve(fx,fy,x,y) %得四个解(0,0),(-1/3,-6),(-7/6,-7/2),(5/6,-5/2).计算Hesse矩阵 fh2=[diff(fx,x),diff(fx,y);diff(fy,x),diff(fy,y)] %计算 eig(subs(fh2,[x,y],[0,0])) %特征值均大于零,得知正定,所以是(0,0)极小值点.极小值用 subs(f,[x,y],[0,0]) %求得9。 %同理可得(-1/3,-6)为极大值点,其它两个为鞍点。 %Exercise 9(以第一小题为例) syms y;f=exp(2*y)/(exp(y)+2); fi=simple(int(f,y)) s=simple(diff(fi)-f) syms a x;f=x^2/sqrt(a^2-x^2); fi=simple(int(f,x)) s=simple(diff(fi)-f) syms a b x;f=1/x/(sqrt(log(x)+a)+sqrt(log(x)+b)); fi=simple(int(f,x)) s=simple(diff(fi)-f) %Exercise 10 syms x y;f=(x-y)^3*sin(x+2*y);Ix=simple(int(f,y,-x,x)) %Exercise 11(1) syms x; f=1/sqrt(2*pi)*exp(-x^2/2); Ix=int(f,x,0,1); vpa(Ix) %Exercise 11(2) syms x; f=exp(2*x)*cos(x)^3; Ix=int(f,x,0,2*pi); vpa(Ix) %Exercise 11(3) syms x; f=x*log(x^4)*asin(1/x); Ix=int(f,x,1,3); vpa(Ix) %Exercise 11(4) syms x; f=sin(x)/x; Ix=int(f,x,0,1); vpa(Ix) %Exercise 11(5) syms r th; f=sqrt(1+r^2*sin(th)); Ix=int(f,r,0,1); Ix=int(Ix,th,0,2*pi); vpa(Ix) %可能死机 %Exercise 11(6) syms x y; f=1+x+y^2; Ix=int(f,y,-sqrt(2*x-x^2),sqrt(2*x-x^2)); Ix=int(Ix,x,0,2); vpa(Ix) %Exercise 12 %1(3) syms x;solve(5*x^23-6*x^7+8*x^6-5*x^2) %1(4) syms x;solve((2*x+3)^3-4) %2 syms x;solve(x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x) %5 syms x y z;[x,y,z]=solve(9*x^2+36*y^2+4*z^2-36,x^2-2*y^2-20*z,16*x-x^3-2*y^2-16*z^2) %6 syms a b;[a,b]=solve(a-0.7*sin(a)-0.2*cos(b),b-0.7*cos(a)+0.2*sin(b)) %7 syms x y;[x,y]=solve((x-2)^2+(y-3+2*x)^2-5,2*(x-3)^2+(y/3)^2-4) vpa(x),vpa(y) %结果虚部接近0, 可舍去 %Exercise 13 %1(1) s=dsolve('Dy=x+y','y(0)=1','x');subs(s,'x',1:3) %1(2) [x,y]=dsolve('Dx=-2*x-3*y','Dy=2*x+y','x(0)=-2.7','y(0)=2.8','t') ezplot(x,y,[0 10]) %1(3) dsolve('D2y-0.01*Dy^2+2*y=sin(t)','y(0)=0','Dy(0)=1','t')(解不出) %1(4) x=dsolve('2*D2x-5*Dx+3*x=45*exp(2*t)','x(0)=2','Dx(0)=1','t') ezplot(x,[0,2]) %1(5) x=dsolve('D2y+(y^2-1)*Dy+y=0','y(0)=2','Dy(0)=0','x')(解不出) %1(6) x=dsolve('D2x=-2/t*Dx+(2*x+10*cos(log(t)))/t/t','x(1)=1','x(3)=3','t'); subs(x,'t',[1.5,2,2.5]) %2 y=dsolve('D5y+10*D4y+54*D3y+132*D2y+137*Dy+50*y','t') %3 [x,y]=dsolve('Dx=-1000.25*x+999.75*y+0.5','Dy=999.75*x-1000.25*y+0.5','x(0)=1','y(0)=-1','t') %Exercise 14 %6(1) ezplot('x^2*sin(x^2-x-2)',[-2 2]) %6(2) ezplot('x^2/4+y^2/9=1') %6(3) ezmesh('x^2+y^2',[-3 3 -3 3]) %6(4) ezmesh('x^4+3*x^2+y^2-2*x-2*y-2*x^2*y+6',[-3 3 -3 13]) %6(5) ezplot3('sin(t)','cos(t)','cos(2*t)',[0 2*pi]) %6(6) ezmesh('2*sin(ph)*cos(th)','2*sin(ph)*sin(th)','2*cos(ph)',[0 pi/2 0 2*pi]) %6(7) ezplot('x','sin(x)',[0 pi]);hold on; ezplot('x','sin(x)*sin(10*x)',[0 pi]); ezplot('x','-sin(x)',[0 pi]); hold off; %Exercise 15 maple('mtaylor((x^2-2*x)*exp(-x^2-y^2-x*y),[x=0,y=a],3)') %Exercise 16(1) x=1:100;p=prod(x) %数值方法 x=sym(1:100);p=prod(x) %符号方法 %Exercise 16(2) x=1:200;p=prod(x) %数值方法 x=sym(1:200);p=prod(x),vpa(p) %符号方法 %Exercise 17 %本题涉及符号和数值混合编程,为了结果简洁,保留5位有效数字。 %先写求Fourier系数的程序 function [a,b]=four(fun,L,n) syms x; a=1/L*int(fun*cos(n*pi*x/L),x,-L,L); b=1/L*int(fun*sin(n*pi*x/L),x,-L,L); a=vpa(a,5);b=vpa(a,5); %下面求解(原题补充说明在[0,2*pi]展开) clear;digits(5); syms x;fun=x*(x-pi)*(x-2*pi) a=0;b=2*pi;L=(b-a)/2; f=1/2*four(fun,L,0); for n=1:6 [a,b]=four(fun,L,i); f=f+a*cos(n*pi*x/L)+b*sin(n*pi*x/L); end f 更多相关推荐: 《数学实验》报告matlab_第六次作业数学实验报告实验名称matlab常微分概率和统计作图学院机械工程学院专业班级姓名学号20xx年11月一实验目的掌握运用Matlab解常微分方程的方法以及运用数值法解常微分的方法和步骤学会使用matlab一组数据... Matlab数学实验报告模板数学实验报告日期20年月日Maylab实验报告学院数学系班级09级A班姓名学号名称二0一二年五月八日 《数学实验》报告matlab-_第三次作业数学实验报告实验名称matlab绘图学院机械工程学院20xx年10月一实验目的掌握Matlab绘图的基本知识学会使用matlab绘制三维曲线和三维曲面掌握基本的绘图指令学会为图形添加各种标注以及改变图形的着色效... matlab数学实验报告2西安交通大学数学实验报告数学实验报告20xx年6月12日2西安交通大学数学实验报告培养容器温度变化率模型一实验目的利用matlab软件估测培养容器温度变化率二实验问题现在大棚技术越来越好能够将温度控制在一定温度... matlab――大学数学实验报告济南大学20xx20xx学年第二学期数学实验上机考试题班级计科1201学号20xx12220xx姓名黄静考试时间20xx年6月17日授课教师王新红说明每题分值20分第5题第6题第7题和第8题可以任选其一第9题和... matlab11实验报告辽宁工程技术大学上机实验报告1程序x2025303540455055606539Xones101xY13215116417117918719621222524339bbintrrintstatsregressY... Matlab实验报告五(微分方程求解Euler折线法)数学与信息科学系实验报告实验名称微分方程求解所属课程数学软件与实验实验类型综合型实验专业信息与计算科学班级学号姓名指导教师1234 东南大学几何与代数matlab实验报告数学实验报告学号姓名吴雪松得分实验1求解线性方程组实验内容用MATLAB求解如下线性方程组Axb其中5100A00006510000006510000006510000006510000006510000006... 数学实验“微分方程组数值算法――四阶Runge-Kutta数值算法”实验报告(内含matlab程序)西京学院数学软件实验任务书实验二十六实验报告一实验名称微分方程组数值算法四阶RungeKutta数值算法二实验目的进一步熟悉微分方程组数值算法四阶RungeKutta数值算法三实验要求运用MatlabCCJav... 数学实验“几种常见的求积分近似解的方法”实验报告(内含matlab程序)西京学院数学软件实验任务书实验二十一实验报告一实验名称Romberg积分法Gauss型积分法高斯勒让德积分法高斯切比雪夫积分法高斯拉盖尔积分法高斯埃尔米特积分法二实验目的进一步熟悉Romberg积分法Gauss... matlab实验报告MATLAB教程及实训班级学号姓名实验报告实验1熟悉matlab的开发环境及矩阵操作一实验的教学目标通过本次实验使学生熟悉MATLAB70的开发环境熟悉MATLAB工作界面的多个常用窗口包括命令窗口历史命令窗口... matlab基本操作实验报告实验一matlab基本操作一实验目的熟悉matlab的安装与启动熟悉matlab用户界面熟悉matlab功能建模元素熟悉matlab优化建模过程二实验设备与仪器1微机2matlab仿真软件三实验步骤1了解mat... matlab数学实验报告(22篇) |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |