MATLAB實(shí)驗(yàn)指導(dǎo)書(shū)加程序+上機(jī)實(shí)例

1、.專(zhuān)業(yè)整理 .MATLAB語(yǔ)言實(shí)驗(yàn)指導(dǎo)書(shū)華東交通大學(xué)電氣學(xué)院永賢2006年2月. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .實(shí)驗(yàn)一MATLAB工作環(huán)境熟悉及簡(jiǎn)單命令的執(zhí)行一、實(shí)驗(yàn)?zāi)康模?熟悉 MATLAB的工作環(huán)境，學(xué)會(huì)使用MATLAB進(jìn)行一些簡(jiǎn)單的運(yùn)算。二、實(shí)驗(yàn)容： MATLAB的啟動(dòng)和退出，熟悉 MATLAB的桌面（ Desktop ），包括菜單（ Menu）、工具條 （ Toolbar ）、命令窗口 (CommandWindow)、歷史命令窗口、 工作空間 (Workspace)等；完成一些基本的矩陣操作；學(xué)習(xí)使用在線幫助系統(tǒng)。三、實(shí)驗(yàn)步驟：1、啟動(dòng) MATLAB，熟悉 MATLAB的桌面。2、在命令窗口

2、執(zhí)行命令完成以下運(yùn)算，觀察workspace 的變化，記錄運(yùn)算結(jié)果。（ 1）（ 365-52 2-70 ） 3（ 2） area=pi*2.52（ 3）已知 x=3， y=4，在 MATLAB中求 z：x2 y3zxy 2（4）將下面的矩陣賦值給變量m1,在 workspace 中察看 m1在存中占用的字節(jié)數(shù)。162313511108m1=97612414151執(zhí)行以下命令m1( 2 , 3 )m1( 11 )m1( : , 3 )m1( 2 : 3 , 1 : 3 )m1( 1 ,4 ) + m1( 2 ,3 ) + m1( 3 ,2 ) + m1( 4 ,1)（ 5）執(zhí)行命令 help a

3、bs查看函數(shù) abs 的用法及用途，計(jì)算abs( 3 + 4i )（ 6）執(zhí)行命令x=0:0.1:6*pi;y=5*sin(x);plot(x,y)（ 6）運(yùn)行 MATLAB的演示程序， demo，以便對(duì) MATLAB有一個(gè)總體了解。四、思考題1、以下變量名是否合法？為什么？（ 1） x2（ 2） 3col（ 3） _row（ 4） for2、求以下變量的值，并在MATLAB中驗(yàn)證。（ 1） a = 1 : 2 : 5 ;（ 2） b = a a a ;（ 3） c = a + b ( 2 , : ). 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .實(shí)驗(yàn)二MATLAB語(yǔ)言矩陣運(yùn)算一、實(shí)驗(yàn)?zāi)康模?掌握基本的矩陣運(yùn)算及

4、常用的函數(shù)。二、實(shí)驗(yàn)容：12324110a56b3c41521、下列運(yùn)算是否合法，為什么？如合法，結(jié)果是多少？(1) result1 = a(2) result2 = a * b(3) result3 = a + b(4) result4 = b * d(5) result5 = b ; c * d(6) result6 = a . * b(7) result7 = a . / b(8) result8 = a . * c(9) result9 = a . b(10) result10 = a . 2(11) result11 = a 2(12) result11 = 2 . a2、用 MAT

5、LAB求下面的的方程組。7212x1491532x27(1)2115x31213213x40xyz1x2yzw8(2)2xy3w33x3y5z6w57212915323、已知 A221151 3213(1) 求矩陣 A 的秩 (rank)(2) 求矩陣 A 的行列式 (determinant)(3) 求矩陣 A 的逆 (inverse)(4) 求矩陣 A 的特征值及特征向量(eigenvalue and eigenvector)4、關(guān)系運(yùn)算與邏輯運(yùn)算已知 a=20,b=-2,c=0,d=1(1) r1 = a b(2) r2 = a b & c d147d 8 5 23 6 0. 學(xué)習(xí)幫手 .

6、專(zhuān)業(yè)整理 .(3) r3 = a = b* (-10)(4) r4 = b | c三、思考題10y2n2 102 929210, 求 y=？（用 format long查看 y 的值）n10實(shí)驗(yàn)三程序的編輯及調(diào)試一、實(shí)驗(yàn)?zāi)康模?掌握 MATLAB程序編輯、運(yùn)行及調(diào)試方法。二、實(shí)驗(yàn)容：1 、 啟 動(dòng)MATLAB 后 ， 點(diǎn) 擊File|New|M-File， 啟 動(dòng)MATLAB 的 程 序 編 輯 及 調(diào) 試 器（ Editor/Debugger ），編輯以下程序， 點(diǎn)擊 File|Save 保存程序， 注意文件名最好用英文字符。點(diǎn)擊 Debug|Run 運(yùn)行程序，在命令窗口查看運(yùn)行結(jié)果，程序如

7、有錯(cuò)誤則改正。注：數(shù)論中一個(gè)有趣的題目：任意一個(gè)正整數(shù)，若為偶數(shù)，則用2除之，若為奇數(shù)，則與3相乘再加上 1。重復(fù)此過(guò)程，最終得到的結(jié)果為1。如：21310516842163105168421運(yùn)行下面的程序，按程序提示輸入n=1,2,3,5,7等數(shù)來(lái)驗(yàn)證這一結(jié)論。%classic 3n+1 problem from number theory.while1n=input(Enter n,negative quits:);ifn1ifrem(n,2)=0n=n/2;elsen=3*n+1;enda=a,n;endaendm2、編程求滿(mǎn)足i 1 2i10000 的最小 m值。三、思考題用對(duì)分法求解

8、方程2e xsin x 在0 ， 1 的解，并驗(yàn)證，在程序中統(tǒng)計(jì)出對(duì)分次數(shù)。提示：先將原方程轉(zhuǎn)化成f ( x) 2e xsin x0 的形式。對(duì)分法的基本思想是：一個(gè)一元方程f(x)=0 ，若 f(x1)*f(x2)0，則在 x1,x2區(qū)間有實(shí)數(shù)解。取該區(qū)間的中點(diǎn)xm=(x1+x2)/2 ，判定 f(x1)和 f(x2) 二者中哪一個(gè)與 f(xm)異號(hào)，. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .若 f(x1)*f(xm)v=myvander(2 3 4 5)得 v=1111234549162582764125生成一些數(shù)據(jù)測(cè)試你寫(xiě)的函數(shù)。三、思考題. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .編寫(xiě)程序，用如下迭代公式求a， a

9、的值分別為：3， 17， 113。迭代的終止條件為x n 1x n10 5，迭代初值x 01.0, 迭代次數(shù)不超過(guò)100 次。分別對(duì)迭代結(jié)果和準(zhǔn)確值進(jìn)行比較，并統(tǒng)計(jì)迭代次數(shù)。x2a2x2x2ax n 1x n 1a22x n 1實(shí)驗(yàn)五 MATLAB的繪圖1、在同一坐標(biāo)系下繪制下面三個(gè)函數(shù)在t 0 ， 4 的圖象。y 1ty 2ty 34e0 . 1 tsin( t )2、編寫(xiě)程序，選擇合適的步距，繪制下面函數(shù)在區(qū)間-6 ， 6 中的圖象。sin x,x0y( x)x,0x 3x 6,x33、用 compass 函數(shù)畫(huà)下面相量圖ua = 1 ; ub = cos(-2*pi/3)+sin(-2*

10、pi/3)*i ; uc=cos(2*pi/3)+sin(2*pi/3)*i;compass(ua,ub,uc,ua-ub,ub-uc,uc-ua)902120601.51501300.518002103302403002704、三維空間曲線繪制z=0:0.1:4*pi;x=cos(z);. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .y=sin(z);plot3(x,y,z)5、用 mesh或 surf函數(shù)，繪制下面方程所表示的三維空間曲面，x 和 y 的取值圍設(shè)為-3 ，3 。x2y2z1010三、思考題在同一坐標(biāo)系下，用不同顏色和線型繪制以下兩個(gè)函數(shù)在t -2 ， 2 圍的圖象。y120.5 ty 22e

11、0.2 t實(shí)驗(yàn)六 MATLAB數(shù)值運(yùn)算一、實(shí)驗(yàn)?zāi)康模?掌握 MATLAB常用的數(shù)值運(yùn)算函數(shù)。二、實(shí)驗(yàn)容：1、求代數(shù)方程 3x54x47x32x29x 120 的 5 個(gè)根，并將其用星號(hào)（* ）標(biāo)記在復(fù)平面圖上。（用 roots和 plot函數(shù)）。2、求代數(shù)方程x510 的 5個(gè)根，并將其用星號(hào)（* ）標(biāo)記在復(fù)平面圖上。 （用 roots和plot函數(shù)）。3、求下面函數(shù)在0.5,4區(qū)間的過(guò)零點(diǎn)。 （用 fzero函）f (x) x32x2 sin(x)5x cos(x)1x4、已知 R=50歐姆， U=4V，二極管D 正向電流與電壓的關(guān)系為：U d q1I dI se KT其中：U 為二極管正向

12、電壓dI s 為反向飽合電流，取 10-12 AK 為玻爾茨曼常數(shù)， 1.38*10 -23T 為絕對(duì)溫度，取300 開(kāi)爾文（ 27 攝氏度）. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .q 為電子電荷 1.6*10 -19 C求此電路中的電流I d 和二極管正向電壓Ud（要求用 fsolve 函數(shù)求解）5、實(shí)驗(yàn)數(shù)據(jù)處理：已知某壓力傳感器的測(cè)試數(shù)據(jù)如下表p0.01.12.12.84.25.06.16.98.19.09.9u1011131417182224293439p 為壓力值， u 為電壓值， 試用多項(xiàng)式 u ( p )ap 3bp 2cpd 來(lái)擬合其特性函數(shù)，求出a,b,c,d，并把擬合曲線和各個(gè)測(cè)試數(shù)據(jù)點(diǎn)

13、畫(huà)在同一幅圖上。實(shí)驗(yàn)七 MATLAB應(yīng)用1、以原點(diǎn)為奇對(duì)稱(chēng)中心的方波y(wt) ，可以用相應(yīng)頻率的基波及其奇次諧波合成。y(wt)4111sin(2n 1)wtsin wtsin 3wtsin 5wt(2n351)n 1,2,3,取的階數(shù)越多， 越接近方波，但總消除不了邊緣上的尖峰，這稱(chēng)為吉布斯效應(yīng)。設(shè)方波頻率為50Hz, 時(shí)間 t 取 00.04 秒 (f=50Hz,w=2*pi*f,h=1e-5,tf=40e-3,t= 0:h:tf) ，編寫(xiě)程序，畫(huà)出如下用 1次諧波、 1,3 次諧波、 1,3,5,7,9次諧波， 1,3,5,19 次諧波合成的近似方波。( 產(chǎn)生方波的函數(shù)為 :squar

14、e)2、用 Simulink求解下圖所示電路 0100 微秒的響應(yīng)。已知 R=6*10-4 歐， C=1700微法， L=6*10-9 享， uc(0)=15kV 。RiGain2Scope1/L1-1/C1UcisisUcULGain1IntegratorGainIntegrator1Scope1模塊參數(shù)設(shè)置 :Integrator1的 Initial condition:15kV在命令窗口為 R,L,C 賦值。仿真參數(shù)設(shè)置如下：Start time:0. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .Stop time:100e-6Solver Type:Variable-stepSolver:ode45Max

15、 step size:1e-7Min step size:autoInitial step size:autoRelative tolerance:1e-3Absolute tolerance:1e-6MATLAB實(shí)驗(yàn)程序?qū)嶒?yàn) 1第1題.(1)x=2*sin(85*pi/180)/(1+exp(2)x =0.2375(2).x=2 1+2i;-0.45 5;y=0.5*log(x+sqrt(1+x(2)y =0.7114 - 0.0253i 0.8968 + 0.3658i0.2139 + 0.9343i 1.1541 - 0.0044i或x=2 1+2i;-0.45 5;d=0.5*log(

16、x+sqrt(1+x*x)d =0.7114 - 0.0253i 0.8968 + 0.3658i0.2139 + 0.9343i 1.1541 - 0.0044i或x=2 1+2*i;-0.45 5;y=0.5*log(x+sqrt(1+x(2)y =0.7114 - 0.0253i 0.8968 + 0.3658i0.2139 + 0.9343i 1.1541 - 0.0044i(3).a=-3.0:0.1:3.0;g=(exp(0.3*a)-exp(-0.3*a).*sin(a+0.3)/2+log(0.3+a)/2)結(jié)果略(4) t=0:0.5:2.5; f1=t.2; f2=t.2-

17、1; f3=t.2-2*t+1; z=(t=0&t=1&t=2&t A=12 34 -4;34 7 87;3 65 7;. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 . B=1 3 -1;2 0 3;3 -2 7; A+6*Bans =18 52 -1046 7 10521 53 49 A-B+eye(3) ans =12 32 -233 8 851681（ 2） A*Bans =68 44 62309 -72 596154 -5 241 A.*Bans =12102468 0 2619 -130 49(3) A3ans =37226 23382448604247370 149188 60076678688 45

18、4142 118820 A.3ans =172839304-643930434365850327274625343(4) A/Bans =16.4000 -13.6000 7.600035.8000 -76.2000 50.200067.0000 -134.0000 68.0000 BAans =109.4000 -131.2000 322.8000-53.000085.0000 -171.0000-61.600089.8000 -.2000(5) A,Bans =1234-413-13478720336573-27 A(1,3,:);B2 ans =1234-436574511101920-

19、540第 3 題. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 . A=1 2 3 4 5;6 7 8 9 10;11 12 13 14 15;16 17 18 19 20;21 22 23 24 25; B=3 0 16;17 -6 9;0 23 -4;9 7 0;4 13 11; C=A*BC =9315077258335237423520397588705557753890717 D=C(3:5,2:3)D =520397705557890717第 4 題(1) a=100:999; k=find(rem(a,21)=0);%找出能杯 21 整除的元素位置 ,find()函數(shù)找出不為0 的元素位置 x=len

20、gth(k)%獲得向量 k 的元素個(gè)數(shù)并賦值給變量 xx =43 k=find(rem(a,21)=0) %顯示能杯 21 整除的元素位置k =Columns 1 through 24627486990111132153174195216237258279300321342363384405426447468489Columns 25 through 43510531552573594615636657678699720741762783804825846867888 y=100+k-1%顯示能杯 21整除的元素y =Columns 1 through 23105126147168189210

21、231252273294315336357378399420441462483504525546567Columns 24 through 43588609630651672693714735756777798819840861882903924945966987(2)sh=CDe345Efg69K; k=find(sh=A&sh sh(k)=;%刪除大寫(xiě)字母 x=sh(1:end)%顯示處理后的字符x =e345fg69實(shí)驗(yàn) 2第 1 題a=1 2 3;4 5 6;b=2 4 -1;1 3 5;c=1;0;-2;d=1 4 7;8 5 2;3 6 0; result1=a%a 的轉(zhuǎn)置resu

22、lt1 =. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .1 42 53 6 result2=a*b%error 應(yīng)采用點(diǎn)乘 result3=a+b%求兩個(gè)矩陣的和result3 =3625811 result4=b*d%矩陣相乘result4 =312222404913 result5=b;c*d result5 =31 22 2240 49 13 -5-87 result6=a.*b%矩陣點(diǎn)乘result6 =28-341530 result7=a./b%矩陣右點(diǎn)除result7 =0.50000.5000-3.00004.00001.66671.2000 result8=a.*c%errora 和 c 維數(shù)

23、不同 result9=a.b%矩陣左點(diǎn)除result9 =2.00002.0000-0.33330.25000.60000.8333 result10=a.2result10=1 4 916 25 36 result11=a2%error等價(jià)于 a*a result12=2.a result12 =2 4 816 32 64第 2 題(1) A=7 2 1 -2;9 15 3 -2;-2 -2 11 5;1 3 2 13; B=4;7;-1;0; X=inv(A)*B%AB等價(jià)于 inv(A)*B,A/B 等價(jià)于 A*inv(B)X =0.49790.14450.0629-0.0813 X1=

24、ABX1 =0.49790.14450.0629 -0.0813(2) a=1 1 1 0;1 2 1 -1;2 -1 0 -3;3 3 5 -6; b=1;8;3;5;. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 . x=inv(a)*bx = 1.00005.0000 -5.0000 -2.0000第 3 題 A=7 2 1 -2;9 15 3 -2;-2 -2 11 5;1 3 2 13; a1=rank(A)a1 =4 a2=det(A)a2 =12568 a3=inv(A)a3 =0.1744 -0.0303 -0.0125 0.0270 -0.1050 0.0789 -0.0121 0.0006 0.

25、0083 0.0173 0.0911 -0.03110.0095 -0.0185 -0.0103 0.0795 V,D=eig(A)%V為向量 A 的特征向量 ,D 為特征值V =-0.76290.0919 + 0.0640i0.0919 - 0.0640i -0.02990.62230.6087 + 0.0276i0.6087 - 0.0276i 0.26370.0807-0.7474-0.74740.6434-0.15540.0342 - 0.2374i0.0342 + 0.2374i 0.7180D =4.8554000012.6460 + 1.8333i000012.6460 - 1.

26、8333i000015.8526第 4 題 a=20; b=-2; c=0; d=1; r1=ab r1 =1 r2=ab&cd r2 =0 r3 = a = b* (-10)r3 =1 r4 = b | c r4 =05 y=0;for k=-10:10 y=y+pow2(k); endformat long yy =2.047999023437500e+003. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .實(shí)驗(yàn) 3方法一a=0;for i=1:20a=a+pow2(i);if a10000m=i;breakendendm方法二a=0;i=1;while (a10000)a=a+pow2(i);i=i+1;en

27、dm=i-1;m實(shí)驗(yàn) 4第 1 題function y=myfun1(x)if x0&x3y=-x+6;end運(yùn)行結(jié)果 : y=myfun1(-pi/2)y = -1 y=myfun1(0)y = 0 y=myfun1(2)y =2 y=myfun1(4)y =2第 2 題function m_x,max_x,min_x,rms_x=myfun2(x) %求平均值sum_x=sum(x);%向量元素求和m,n=size(x);%最好用 n=length(x);m_x=sum_x/n;%求最大值采用逐個(gè)比較方式if x(1)x(2)max_x=x(1);elsemax_x=x(2);endfor

28、 k=3:nif max_xx(k). 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .max_x=x(k);elsemax_x=max_x;%可省略endend%求最小值if x(1)x(k)min_x=x(k);elsemin_x=min_x;%可省略endend%求均方根值sum_x2=0;for k=1:nsum_x2=sum_x2+x(k).2;rms_x=sqrt(sum_x2/n);endm_x;max_x;min_x;rms_x;%按照函數(shù)值行參順序輸出結(jié)果運(yùn)行結(jié)果 : m_x,max_x,min_x,rms_x=myfun2(sin(0:0.01:6*pi) m_x =-1.1256e-007 ma

29、x_x =1.0000 min_x =-1.0000 rms_x =0.7071 m_x,max_x,min_x,rms_x=myfun2(rand(1,200) m_x =0.4977 max_x =0.9961min_x =0.0046rms_x =0.5778第 3 題function v=myvander(x)v1=vander(x);%生成德蒙矩陣v2=v1;v=flipud(v2);%實(shí)現(xiàn)矩陣上下翻轉(zhuǎn)運(yùn)行結(jié)果 : v=myvander(2:5)v =. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .111123454 916258 27 64 125思考題function x,n=sqrt_a(a)x=

30、1.0;for k=1:100m=x;x=x/2+a/(2*x);if abs(x-m)=10(-5)breakendendx;n=k;s=(x-sqrt(a);if s x,n=sqrt_a(3)正確x = 1.7321n = 5 x,n=sqrt_a(17)正確x = 4.1231n = 6 x,n=sqrt_a(113)正確x = 10.6301n = 8實(shí)驗(yàn) 5第1題.方法 1 t=linspace(0,4*pi,200); y1=t;y2=sqrt(t); y3=4*pi*exp(-0.1*t).*sin(t);plot(t,y1,b,t,y2,g,t,y3,r)方法 2 t=lin

31、space(0,4*pi,200);y1=t;y2=sqrt(t); y3=4*pi*exp(-0.1*t).*sin(t); t=t,t,t;y=y1,y2,y3;plot(t,y). 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .第 2 題 x=linspace(-6,6,100); y=;for x0=xif x0=0 y=y,sin(x0);elseif x0 z=0:0.1:4*pi;x=cos(z);y=sin(z);plot3(x,y,z,rp);title( 三維空間曲線 ); text(0,0,0,origin); xlabel(X),ylabel(Y),zlabel(Z); grid;第 5 題

32、（1） x=-3:0.1:3; x,y=meshgrid(x); z=-x.2/10+y.2/10; mesh(x,y,z); xlabel(X),ylabel(Y),zlabel(Z); title( 立體網(wǎng)狀圖 );（2） x=-3:0.1:3; x,y=meshgrid(x); z=-x.2/10+y.2/10; surf(x,y,z);xlabel(X),ylabel(Y),zlabel(Z);title(立體曲面圖 );思考題 t=-2*pi:pi/100:2*pi; y1=2.(0.5*abs(t); y2=2*exp(-0.2*t); plot(t,y1,-g);hold on;

33、 plot(t,y2,:r);legend( 曲線 y1, 曲線 y2); hold off;grid;. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .實(shí)驗(yàn) 6第 1 題 A=3,4,7,2,9,12;x=roots(A)plot(x,*);grid; x =-0.8612 + 1.4377i -0.8612 - 1.4377i 0.6737 + 1.0159i0.6737 - 1.0159i -0.9583第 2 題 A=1,0,0,0,0,-1; x=roots(A) plot(x,*);grid; x =-0.8090 + 0.5878i -0.8090 - 0.5878i 0.3090 + 0.9511i

34、0.3090 - 0.9511i1.0000第 3 題%估計(jì)零點(diǎn) fplot(x3+1/x,0.5,4); hold on; fplot(2*x2*sin(x)-5*x*cos(x),0.5,4); hold off;m,n=ginput%建立函數(shù)function y=f(x)y=x3-2*x2*sin(x)+5*x*cos(x)+1/x;%調(diào)用函數(shù) y1=fzero(fz,1.5)y1 =1.5117 y2=fzero(fz,2.5)y2 =2.6095第 4 題%估計(jì)零點(diǎn) axis(0,1,0,0.1);fplot(10(-12)*exp(Ud*1.6*10(-19)/(1.38*10(-

35、23)*300)-1),0,4); hold on;fplot(Ud-4)/50,0,4); hold off; m,n=ginput%建立函數(shù)function f=myfun(X)Id=X(1);Ud=X(2);f(1)=Id-10(-12)*exp(Ud*1.6*10(-19)/(1.38*10(-23)*300)-1);f(2)=50*Id-Ud-4;. 學(xué)習(xí)幫手 .專(zhuān)業(yè)整理 .%調(diào)用函數(shù) x=fsolve(myfun,0,0.8,optimset(Display,off)x =0.09360.6795%驗(yàn)證結(jié)果 K=myfun(x) K =1.0e-008 *-0.30870第 5 題

36、 p=0.0,1.1,2.1,2.8,4.2,5.0,6.1,6.9,8.1,9.0,9.9;u=10,11,13,14,17,18,22,24,29,34,39;x=polyfit(p,u,3) %得多項(xiàng)式系數(shù) t=linspace(0,10,100);y=polyval(x,t);%求多項(xiàng)式得值plot(p,u,*,t,y,r)%畫(huà)擬和曲線x =0.0195-0.04121.44699.8267實(shí)驗(yàn) 7第 1 題 f=50; w=2*pi*f; h=1e-5; tf=40e-3; t=0:h:tf; wt=w*t; y1=4/pi*sin(wt); y2=4/pi*(sin(wt)+1/3*sin(3*wt); y3=4/pi*(sin(wt)+1/3*sin(3*wt)+1/5*sin(5*wt)+1/7*sin(7*wt)+1/9*sin(9*wt); y4=4/pi*(sin(wt)+1/3*sin(3*wt)+1/5*sin(5*wt)+1/7*sin(7*wt)+1/9*sin(9*wt)+1/11*s in(11*wt)+1/13*sin(13*wt)+1/15*sin(15*wt)+1/17*sin(17*wt)+1/19*sin(19*w