chap8.bas

介绍一些matlab在数学方面的应用

Attribute VB_Name = "chap8"
Public Sub PutOnClipboard(str As Variant)
    Dim MyDataObj As New DataObject
    MyDataObj.SetText Format(str)
    MyDataObj.PutInClipboard
End Sub
Sub h1()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=0:.12:1;y=(x.^2-3*x+5).*exp(-5*x).*sin(x); " + CR + "plot(x,y,x,y,'o')")
End Sub
Sub h2()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x1=0:.02:1; y0=(x1.^2-3*x1+5).*exp(-5*x1).*sin(x1);" + CR + "y1=interp1(x,y,x1); y2=interp1(x,y,x1,'cubic');" + CR + "y3=interp1(x,y,x1,'spline'); y4=interp1(x,y,x1,'nearest');" + CR + "plot(x1,[y1',y2',y3',y4'],':',x,y,'o',x1,y0)")
End Sub
Sub h3()
PutOnClipboard ("[max(abs(y0(1:49)-y2(1:49))),max(abs(y0-y3)),max(abs(y0-y4))]")
End Sub
Sub h4()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=-1+2*[0:10]/10; y0=1./(1+25*x0.^2);" + CR + "x=-1:.01:1; y=lagrange(x0,y0,x); % Lagrange 插值" + CR + "ya=1./(1+25*x.^2); plot(x,ya,x,y,':')")
End Sub
Sub h5()
CR = Chr(13) + Chr(10)
PutOnClipboard ("y1=interp1(x0,y0,x,'cubic'); " + CR + "y2=interp1(x0,y0,x,'spline');" + CR + "plot(x,ya,x,y1,':',x,y2,'--')")
End Sub
Sub h6()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=0:pi/30:pi; y0=sin(x0); " + CR + "I=quadspln(x0,y0,0,pi)")
End Sub
Sub h7()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=0:pi/10:pi; y0=sin(x0); I1=trapz(x0,y0)")
End Sub
Sub h8()
CR = Chr(13) + Chr(10)
PutOnClipboard ("I=quadspln(x0,y0,0,pi)")
End Sub
Sub h9()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=[0,0.4,1 2,pi]; y0=sin(x0); % 生成样本点" + CR + "plot(x0,y0,x0,y0,'o')" + CR + "I=quadspln(x0,y0,0,pi) % 大约有1% 的相对误差，应该说是相当精确的")
End Sub
Sub h10()
CR = Chr(13) + Chr(10)
PutOnClipboard ("I1=trapz(x0,y0) % 用trapz() 函数将得出很大的相对误差(7.9%)")
End Sub
Sub h11()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=[0:0.01:pi,pi]; y0a=sin(x); " + CR + "y=interp1(x0,y0,x,'spline');" + CR + "plot(x0,y0,x,y,':',x,y0a,x0,y0,'o')")
End Sub
Sub h12()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=[0:3*pi/2/200:3*pi/2]; y=cos(15*x); I=quadspln(x,y,0,3*pi/2)")
End Sub
Sub h13()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=[0:3*pi/2/1000:3*pi/2]; y0=cos(15*x0);" + CR + "y1=interp1(x,y,x0,'spline');" + CR + " plot(x0,y,x0,y1,':')")
End Sub
Sub h14()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[x,y]=meshgrid(-3:.6:3, -2:.4:2);" + CR + "  z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);" + CR + "surf(x,y,z), axis([-3,3,-2,2,-0.7,1.5])")
End Sub
Sub h15()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[x1,y1]=meshgrid(-3:.2:3, -2:.2:2);" + CR + " z1=interp2(x,y,z,x1,y1);" + CR + "  surf(x1,y1,z1), axis([-3,3,-2,2,-0.7,1.5])")
End Sub
Sub h16()
CR = Chr(13) + Chr(10)
PutOnClipboard ("z1=interp2(x,y,z,x1,y1,'cubic'); " + CR + "z2=interp2(x,y,z,x1,y1,'spline');" + CR + "surf(x1,y1,z1), axis([-3,3,-2,2,-0.7,1.5])" + CR + "figure; " + CR + "surf(x1,y1,z2), axis([-3,3,-2,2,-0.7,1.5])")
End Sub
Sub h17()
CR = Chr(13) + Chr(10)
PutOnClipboard ("z=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1); % 新网格各点的函数值" + CR + "surf(x1,y1,abs(z-z1)), axis([-3,3,-2,2,0,0.08])" + CR + "figure; " + CR + "surf(x1,y1,abs(z-z2)), axis([-3,3,-2,2,0,0.025])")
End Sub
Sub h18()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=-3+6*rand(200,1); y=-2+4*rand(200,1);" + CR + "z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); % 生成已知数据" + CR + "plot(x,y,'x') % 样本点的二维分布" + CR + "figure, plot3(x,y,z,'x'), axis([-3,3,-2,2,-0.7,1.5])")
End Sub
Sub h19()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[x1,y1]=meshgrid(-3:.2:3, -2:.2:2);" + CR + "z1=griddata(x,y,z,x1,y1,'cubic'); " + CR + "surf(x1,y1,z1), axis([-3,3,-2,2,-0.7,1.5])" + CR + "z2=griddata(x,y,z,x1,y1,'v4');" + CR + "figure;" + CR + " surf(x1,y1,z2), axis([-3,3,-2,2,-0.7,1.5])")
End Sub
Sub h20()
CR = Chr(13) + Chr(10)
PutOnClipboard ("z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1); % 新网格各点的函数值" + CR + "surf(x1,y1,abs(z0-z1)); axis([-3,3,-2,2,0,0.15])" + CR + "figure;" + CR + " surf(x1,y1,abs(z0-z2)); axis([-3,3,-2,2,0,0.15])")
End Sub
Sub h21()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=-3+6*rand(200,1); y=-2+4*rand(200,1); % 重新生成样本点" + CR + "z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);" + CR + "ii=find((x+1).^2+(y+0.5).^2>0.5^2); % 找出不满足条件的点坐标" + CR + "x=x(ii); y=y(ii); z=z(ii); plot(x,y,'x')" + CR + "t=[0:.1:2*pi,2*pi]; x0=-1+0.5*cos(t); y0=-0.5+0.5*sin(t);" + CR + "line(x0,y0) ")
End Sub
Sub h22()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[x1,y1]=meshgrid(-3:.2:3, -2:.2:2);" + CR + "z1=griddata(x,y,z,x1,y1,'v4');" + CR + "surf(x1,y1,z1), axis([-3,3,-2,2,-0.7,1.5])")
End Sub
Sub h23()
CR = Chr(13) + Chr(10)
PutOnClipboard ("z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);" + CR + "surf(x1,y1,abs(z0-z1)), axis([-3,3,-2,2,0,0.1])")
End Sub
Sub h24()
CR = Chr(13) + Chr(10)
PutOnClipboard ("contour(x1,y1,abs(z0-z1),30);" + CR + " hold on, plot(x,y,'x'); line(x0,y0)")
End Sub
Sub h25()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[x,y,z]=meshgrid(-1:.2:1);" + CR + "V=exp(x.^2.*z+y.^2.*x+z.^2.*y).*cos(x.^2.*y.*z+z.^2.*y.*x);" + CR + "[x0,y0,z0]=meshgrid(-1:0.05:1);" + CR + "V0=exp(x0.^2.*z0+y0.^2.*x0+z0.^2.*y0).*cos(x0.^2.*y0.*z0+z0.^2.*y0.*x0);" + CR + "V1=interp3(x,y,z,V,x0,y0,z0,'spline');err=V1-V0;max(err(:))")
End Sub
Sub h26()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=[0,0.4,1 2,pi]; y0=sin(x0);" + CR + "sp=csapi(x0,y0), fnplt(sp,':'); hold on," + CR + "ezplot('sin(t)',[0,pi]); plot(x0,y0,'o')")
End Sub
Sub h27()
PutOnClipboard ("sp.coefs")
End Sub
Sub h28()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x);" + CR + "sp=csapi(x,y); fnplt(sp)" + CR + "c=[sp.breaks(1:4)' sp.breaks(2:5)' sp.coefs(1:4,:)," + CR + "sp.breaks(5:8)' sp.breaks(6:9)' sp.coefs(5:8,:) ];")
End Sub
Sub h29()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=-3:.6:3; y0=-2:.4:2; [x,y]=ndgrid(x0,y0); % 注意这里只能用ndgrid" + CR + "z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); % 否则生成的z 矩阵顺序有问题" + CR + "sp=csapi({x0,y0},z); fnplt(sp);")
End Sub
Sub h30()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=[0,0.4,1 2,pi]; y0=sin(x0); ezplot('sin(t)',[0,pi]); hold on" + CR + "sp1=csapi(x0,y0); fnplt(sp1,'--'); % 三次分段多项式样条插值" + CR + "sp2=spapi(5,x0,y0); fnplt(sp2,':') % 5 次B 样条插值")
End Sub
Sub h31()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x);" + CR + "ezplot('(x^2-3*x+5)*exp(-5*x)*sin(x)',[0,1]), hold on" + CR + "sp1=csapi(x,y); fnplt(sp1,'--'); sp2=spapi(5,x,y); fnplt(sp2,':')")
End Sub
Sub h32()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms x; f=(x^2-3*x+5)*exp(-5*x)*sin(x);" + CR + " ezplot(diff(f),[0,1]), hold on" + CR + "x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x);" + CR + "sp1=csapi(x,y); dsp1=fnder(sp1,1); fnplt(dsp1,'--')" + CR + "sp2=spapi(5,x,y); dsp2=fnder(sp2,1); fnplt(dsp2,':'); axis([0,1,-0.8,5])")
End Sub
Sub h33()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=-3:.3:3; y0=-2:.2:2; [x,y]=ndgrid(x0,y0);" + CR + "z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);" + CR + "sp=spapi({5,5},{x0,y0},z); dspxy=fnder(sp,[1,1]); fnplt(dspxy)")
End Sub
Sub h34()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms x y; z=(x^2-2*x)*exp(-x^2-y^2-x*y);" + CR + "ezsurf(diff(diff(z,x),y),[-3 3],[-2 2])")
End Sub
Sub h35()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=[0,0.4,1 2,pi]; y=sin(x);" + CR + "sp1=csapi(x,y); a=fnint(sp1,1);" + CR + " xx=fnval(a,[0,pi]); xx(2)-xx(1)")
End Sub
Sub h36()
CR = Chr(13) + Chr(10)
PutOnClipboard ("sp2=spapi(5,x,y); b=fnint(sp2,1);" + CR + " xx=fnval(b,[0,pi]); xx(2)-xx(1)")
End Sub
Sub h37()
CR = Chr(13) + Chr(10)
PutOnClipboard ("ezplot('-cos(t)+2',[0,pi]); hold on" + CR + "fnplt(a,'--'); fnplt(b,':')")
End Sub
Sub h38()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=0:.1:1; y0=(x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);" + CR + "p3=polyfit(x0,y0,3); vpa(poly2sym(p3),10)")
End Sub
Sub h39()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=0:.01:1; ya=(x.^2-3*x+5).*exp(-5*x).*sin(x);" + CR + "y1=polyval(p3,x); plot(x,y1,x,ya,x0,y0,'o')")
End Sub
Sub h40()
CR = Chr(13) + Chr(10)
PutOnClipboard ("p4=polyfit(x0,y0,4); y2=polyval(p4,x);" + CR + "p5=polyfit(x0,y0,5); y3=polyval(p5,x);" + CR + "p8=polyfit(x0,y0,8); y4=polyval(p8,x);" + CR + "plot(x,ya,x0,y0,'o',x,y2,x,y3,x,y4)")
End Sub
Sub h41()
PutOnClipboard ("vpa(poly2sym(p8),5)")
End Sub
Sub h42()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms x; y=(x^2-3*x+5)*exp(-5*x)*sin(x);" + CR + "vpa(taylor(y,9),5)")
End Sub
Sub h43()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=-1+2*[0:10]/10; y0=1./(1+25*x0.^2);" + CR + "x=-1:.01:1; ya=1./(1+25*x.^2);" + CR + "p3=polyfit(x0,y0,3); y1=polyval(p3,x);" + CR + "p5=polyfit(x0,y0,5); y2=polyval(p5,x);" + CR + "p8=polyfit(x0,y0,8); y3=polyval(p8,x);" + CR + "p10=polyfit(x0,y0,10); y4=polyval(p10,x);")
End Sub
Sub h44()
PutOnClipboard ("plot(x,ya,x,y1,x,y2,'-.',x,y3,'--',x,y4,':')")
End Sub
Sub h45()
PutOnClipboard ("syms x; y=1/(1+25*x^2); p=taylor(y,x,10)")
End Sub
Sub h46()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x1=-1:0.01:1; ya=1./(1+25*x1.^2); y1=subs(p,x,x1);" + CR + " plot(x1,y1)")
End Sub
Sub h47()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms x; c=taylor(exp(-2*x),10); c=sym2poly(c); c=c(end:-1:1); x=0:0.01:8;" + CR + "nd=[3:7]; xx=[0,2,2+eps,8]; yy=[0,0,1,1]; plot(xx,yy); hold on" + CR + "for i=1:length(nd)" + CR + "[n,d]=padefcn(c,0,nd(i)); y=polyval(n,x)./polyval(d,x); plot(x,y)" + CR + "end")
End Sub
Sub h48()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=[0,0.2,0.4,0.7,0.9,0.92,0.99,1.2,1.4,1.48,1.5]';" + CR + "y=[2.88;2.2576;1.9683;1.9258;2.0862;2.109;2.1979;2.5409;2.9627;3.155;3.2052];" + CR + "A=[ones(size(x)) exp(-3*x), cos(-2*x).*exp(-4*x) x.^2];" + CR + "c=A\y; c1=c'")
End Sub
Sub h49()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=[0:0.01:1.5]';" + CR + "A1=[ones(size(x0)) exp(-3*x0), cos(-2*x0).*exp(-4*x0) x0.^2];" + CR + "y1=A1*c; plot(x0,y1,x,y,'x')")
End Sub
Sub h50()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=[1.1052,1.2214,1.3499,1.4918,1.6487,1.8221,2.0138,2.2255,2.4596,2.7183,3.6693];" + CR + "y=[0.6795,0.6006,0.5309,0.4693,0.4148,0.3666,0.3241,0.2864,0.2532,0.2238,0.1546];" + CR + "plot(x,y,x,y,'*')")
End Sub
Sub h51()
PutOnClipboard ("x1=log(x); y1=log(y); plot(x1,y1,x1,y1,'*')")
End Sub
Sub h52()
PutOnClipboard ("A=[x1' ones(size(x1'))]; c=[A\y1']'")
End Sub
Sub h53()
PutOnClipboard ("exp(c(2))")
End Sub
Sub h54()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=[0:.1:2]'; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); n=7; A=[];" + CR + "for i=1:n+1, A(:,i)=x.^(n+1-i); end" + CR + "c=A\y; vpa(poly2sym(c),5)")
End Sub
Sub h55()
PutOnClipboard ("x=0:.1:10; y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);")
End Sub
Sub h56()
PutOnClipboard ("f=inline('a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x)','a','x');")
End Sub
Sub h57()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y); xx',res")
End Sub
Sub h58()
CR = Chr(13) + Chr(10)
PutOnClipboard ("ff=optimset; ff.TolFun=1e-20; ff.TolX=1e-15; % 修改精度限制" + CR + "[xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y,[],[],ff); xx',res")
End Sub
Sub h59()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x1=0:0.01:10; y1=f(xx,x1); plot(x1,y1,x,y,'o')")
End Sub
Sub h60()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=0.1:0.1:1;" + CR + "y=[2.3201,2.6470,2.9707,3.2885,3.6008,3.9090,4.2147,4.5191,4.8232,5.1275];")
End Sub
Sub h61()
PutOnClipboard ("a=lsqcurvefit('c8f3',[1;2;2;3],x,y); a'")
End Sub
Sub h62()
PutOnClipboard ("y1=c8f3(a,x); plot(x,y,x,y1,'o')")
End Sub
Sub h63()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms A1 A2 w1 w2 t1 t2 t tau T;" + CR + " x=A1*cos(w1*t+t1)+A2*cos(w2*t+t2);" + CR + "Rxx=simple(limit(int(x*subs(x,t,t+tau),t,0,T)/T,T,inf))")
End Sub
Sub h64()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=0:0.01:5; y1=x.*exp(-4*x).*sin(3*x); " + CR + "y2=x.*exp(-4*x).*cos(3*x);" + CR + "R=corrcoef(y1,y2)")
End Sub
Sub h65()
CR = Chr(13) + Chr(10)
PutOnClipboard ("t=0:0.01:5; x=t.*exp(-4*t).*sin(3*t); y=t.*exp(-4*t).*cos(3*t);" + CR + "N=150; c1=xcorr(x,N); x1=[-N:N]; stem(x1,c1) % 实际绘图点选得更稀疏" + CR + "figure; c1=xcorr(x,y,N); x1=[-N:N]; stem(x1,c1)")
End Sub
Sub h66()
CR = Chr(13) + Chr(10)
PutOnClipboard ("h=0.01; t=0:h:1;" + CR + " x=12*sin(2*pi*t+pi/4)+5*cos(2*pi*4*t); X=fft(x);" + CR + "f=t/h; " + CR + "plot(f(1:floor(length(f)/2)),abs(X(1:floor(length(f)/2))))")
End Sub
Sub h67()
CR = Chr(13) + Chr(10)
PutOnClipboard ("ix=real(ifft(X));" + CR + " plot(t,x,t,ix,':'); norm(x-ix)")
End Sub
Sub h68()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=0:0.002:2; y=exp(-x).*sin(5*x);" + CR + " r=0.05*randn(size(x)); y1=y+r;" + CR + " plot(x,y1)")
End Sub
Sub h69()
CR = Chr(13) + Chr(10)
PutOnClipboard ("b=1.2296e-6*conv([1 4 6 4 1],[1 3 3 1]);" + CR + " a=conv([1,-0.7265],conv([1,-1.488,0.5644],conv([1,-1.595,0.6769],[1,-1.78,0.8713])));" + CR + "x=0:0.002:2;y=exp(-x).*sin(5*x);r=0.05*randn(size(x));y1=y+r;")
End Sub
Sub h70()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[h,w]=freqz(b,a,100); plot(w,abs(h)) % 放大倍数绘制" + CR + "figure; y2=filter(b,a,y1); plot(x,y1,x,y2) % 滤波效果")
End Sub
Sub h71()
CR = Chr(13) + Chr(10)
PutOnClipboard ("f1=figure; f2=figure;" + CR + "for n=5:2:20" + CR + "figure(f1); [b,a]=butter(n,0.1); y2=filter(b,a,y1);" + CR + " plot(x,y2); hold on" + CR + "figure(f2); [h,w]=freqz(b,a,100); plot(w,abs(h)); hold on" + CR + "end")
End Sub
Sub h72()
CR = Chr(13) + Chr(10)
PutOnClipboard ("for wn=0.1:0.1:0.7" + CR + "figure(f1); [b,a]=butter(7,wn); y2=filter(b,a,y1);" + CR + " plot(x,y2); hold on" + CR + "figure(f2); [h,w]=freqz(b,a,100); plot(w,abs(h)); hold on" + CR + "end")
End Sub
Sub h73()
CR = Chr(13) + Chr(10)
PutOnClipboard ("maple('with(numtheory):');" + CR + "f = maple(['cfe:=cfrac(pi,20)'])")
End Sub
Sub h74()
CR = Chr(13) + Chr(10)
PutOnClipboard ("n=maple('nthnumer','cfe',4);" + CR + "d=maple('nthdenom','cfe',4);" + CR + "[vpa(n),vpa(d)]")
End Sub
Sub h75()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms x; fun='sin(x)*exp(-x)/(x+1)^3';" + CR + "maple('with(numtheory):');" + CR + "f = maple(['cfe:=cfrac(' fun ',x,10)'])")
End Sub
Sub h76()
CR = Chr(13) + Chr(10)
PutOnClipboard ("n=collect(maple('nthnumer','cfe',8),x);" + CR + "d=collect(maple('nthdenom','cfe',8),x);" + CR + "[n,d]=numden(n/d); G=n/d; latex(G)" + CR + "G1=maple('nthnumer','cfe',10),x)/maple('nthdenom','cfe',10),x);")
End Sub