chap7.bas

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

Attribute VB_Name = "模块1"
Public Sub PutOnClipboard(str As Variant)
    Dim MyDataObj As New DataObject
    MyDataObj.SetText Format(str)
    MyDataObj.PutInClipboard
End Sub
Sub g1()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms t; u=exp(-5*t)*cos(2*t+1)+5;" + CR + "uu=5*diff(u,t,2)+4*diff(u,t)+2*u")
End Sub
Sub g2()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms t y;" + CR + "y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=','87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'])")
End Sub
Sub g3()
CR = Chr(13) + Chr(10)
PutOnClipboard ("y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=','87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'],'y(0)=3','Dy(0)=2','D2y(0)=0','D3y(0)=0')")
End Sub
Sub g4()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[n,d]=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)));" + CR + " [n,d]")
End Sub
Sub g5()
CR = Chr(13) + Chr(10)
PutOnClipboard ("vpa(-445/26*cos(sym(1))-51/13*sin(1)-69/2+8704/185)")
End Sub
Sub g6()
CR = Chr(13) + Chr(10)
PutOnClipboard ("y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=','87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'],'y(0)=1/2','Dy(pi)=1','D2y(2*pi)=0','Dy(2*pi)=1/5')")
End Sub
Sub g7()
PutOnClipboard ("vpa(y,10)")
End Sub
Sub g8()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms t y; u=sin(t); uu=3*diff(u)+3*u")
End Sub
Sub g9()
CR = Chr(13) + Chr(10)
PutOnClipboard ("y=dsolve('D5y+5*D4y+12*D3y+16*D2y+12*Dy+4*y=3*cos(t)+3*sin(t)','y(0)=0','Dy(0)=0','D2y(0)=0','D3y(0)=0','D4y(0)=0')")
End Sub
Sub g10()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[x,y]=dsolve('D2x+2*Dx=x+2*y-exp(-t)','Dy=4*x+3*y+4*exp(-t)')")
End Sub
Sub g11()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms t; u=exp(-5*t)*cos(2*t+1)+5; % 定义输入函数" + CR + " uu=laplace(5*diff(u,2)+4*diff(u)+2*u); % 对等号右侧整个式子进行变换" + CR + "uu=collect(simple(uu)); latex(uu) % 化简")
End Sub
Sub g12()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms s; G=uu/(s^4+10*s^3+35*s^2+50*s+24);" + CR + "Y=residue(G,s); y=ilaplace(Y); latex(simple(y))")
End Sub
Sub g13()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms t; u=exp(-5*t)*cos(2*t+1)+5; % 定义输入函数" + CR + "uu=laplace(5*diff(u,2)+4*diff(u)+2*u) % 求出并显示等号右侧式子")
End Sub
Sub g14()
CR = Chr(13) + Chr(10)
PutOnClipboard ("y1=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=','87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'],'y(0)=0','Dy(0)=0','D2y(0)=0','D3y(0)=0');" + CR + "simple(y-y1) % 验证两种方法的结果完全一致")
End Sub
Sub g15()
PutOnClipboard ("ezplot(y,[0,10]); axis([0,10,0,0.6])")
End Sub
Sub g16()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms t x" + CR + "x=dsolve('Dx=x*(1-x^2)')")
End Sub
Sub g17()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms t x; " + CR + "x=dsolve('Dx=x*(1-x^2)+1')")
End Sub
Sub g18()
PutOnClipboard ("syms y mu; y=dsolve('D2y+mu*(y^2-1)*Dy+y')")
End Sub
Sub g19()
CR = Chr(13) + Chr(10)
PutOnClipboard ("t_final=100; x0=[0;0;1e-10];" + CR + "[t,x]=ode45('lorenzeq',[0,t_final],x0);plot(t,x)," + CR + "figure; % 打开新图形窗口" + CR + "plot3(x(:,1),x(:,2),x(:,3));" + CR + "axis([10 42 -20 20 -20 25]); % 根据实际数值手动设置坐标系")
End Sub
Sub g20()
CR = Chr(13) + Chr(10)
PutOnClipboard ("f1=inline(['[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);','-x(1)*x(2)+28*x(2)-x(3)]'],'t','x');")
End Sub
Sub g21()
CR = Chr(13) + Chr(10)
PutOnClipboard ("t_final=100; x0=[0;0;1e-10]; [t,x]=ode45(f1,[0,t_final],x0);" + CR + "plot(t,x), figure;" + CR + "plot3(x(:,1),x(:,2),x(:,3)); axis([10 42 -20 20 -20 25]);")
End Sub
Sub g22()
CR = Chr(13) + Chr(10)
PutOnClipboard ("t_final=100; x0=[0;0;1e-10];" + CR + "b1=8/3; r1=10; s1=28; % 并不要求变量名必须是beta 等" + CR + "[t,x]=ode45('lorenz1',[0,t_final],x0,[],b1,r1,s1);" + CR + "plot(t,x)," + CR + "figure; plot3(x(:,1),x(:,2),x(:,3)); axis([10 42 -20 20 -20 25]);")
End Sub
Sub g23()
CR = Chr(13) + Chr(10)
PutOnClipboard ("t_final=100; x0=[0;0;1e-10];" + CR + "b2=2; r2=5; s2=20;" + CR + "[t2,x2]=ode45('lorenz1',[0,t_final],x0,[],b2,r2,s2);" + CR + "plot(t2,x2)," + CR + "figure; plot3(x2(:,1),x2(:,2),x2(:,3)); axis([0 72 -20 22 -35 40]);")
End Sub
Sub g24()
CR = Chr(13) + Chr(10)
PutOnClipboard ("f2=inline(['[-beta*x(1)+x(2)*x(3); -rho*x(2)+rho*x(3);','-x(1)*x(2)+sigma*x(2)-x(3)]'],'t','x','flag','beta','rho','sigma');")
End Sub
Sub g25()
CR = Chr(13) + Chr(10)
PutOnClipboard ("t_final=100; x0=[0;0;1e-10];" + CR + "b2=2; r2=5; s2=20;" + CR + "[t2,x2]=ode45(f2,[0,t_final],x0,[],b2,r2,s2);" + CR + "plot(t2,x2)," + CR + "figure; plot3(x2(:,1),x2(:,2),x2(:,3)); axis([0 72 -20 22 -35 40]);")
End Sub
Sub g26()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=[-0.2; -0.7]; t_final=20;" + CR + "mu=1; [t1,y1]=ode45('vdp_eq',[0,t_final],x0,[],mu);" + CR + "mu=2; [t2,y2]=ode45('vdp_eq',[0,t_final],x0,[],mu);" + CR + "plot(t1,y1,t2,y2,':')" + CR + "figure; plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),':')")
End Sub
Sub g27()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=[2;0]; t_final=3000;" + CR + "mu=1000; [t,y]=ode45('vdp_eq',[0,t_final],x0,[],mu);")
End Sub
Sub g28()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=[1.2; 0; 0; -1.04935751];" + CR + "tic, [t,y]=ode45('apolloeq',[0,20],x0); toc" + CR + "length(t), plot(y(:,1),y(:,3))")
End Sub
Sub g29()
CR = Chr(13) + Chr(10)
PutOnClipboard ("options=odeset; options.RelTol=1e-6;" + CR + "tic, [t1,y1]=ode45('apolloeq',[0,20],x0,options); toc" + CR + "length(t1), plot(y1(:,1),y1(:,3)),")
End Sub
Sub g30()
CR = Chr(13) + Chr(10)
PutOnClipboard ("min(diff(t1))")
End Sub
Sub g31()
CR = Chr(13) + Chr(10)
PutOnClipboard ("plot(t1(1:end-1),diff(t1))")
End Sub
Sub g32()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=[1.2; 0; 0; -1.04935751];" + CR + "tic, [t1,y1]=rk_4('apolloeq',[0,20,0.01],x0); toc" + CR + "plot(y1(:,1),y1(:,3)) % 绘制出轨迹曲线")
End Sub
Sub g33()
CR = Chr(13) + Chr(10)
PutOnClipboard ("tic, [t2,y2]=rk_4('apolloeq',[0,20,0.001],x0); toc" + CR + "plot(y2(:,1),y2(:,3)) % 绘制出轨迹曲线")
End Sub
Sub g34()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms x1 x2 x3 x4" + CR + "[dx,dy]=solve('dx+2*x4*x1=2*dy','dx*x4+3*x2*dy+x1*x4-x3=5','dx,dy');" + CR + "latex(dx), latex(dy)")
End Sub
Sub g35()
CR = Chr(13) + Chr(10)
PutOnClipboard ("h_opt=odeset; h_opt.RelTol=1e-6; x0=[2;0]; t_final=3000;" + CR + "tic, mu=1000; [t,y]=ode15s('vdp_eq',[0,t_final],x0,h_opt,mu); toc")
End Sub
Sub g36()
PutOnClipboard ("plot(t,y(:,1)); figure; plot(t,y(:,2))")
End Sub
Sub g37()
CR = Chr(13) + Chr(10)
PutOnClipboard ("syms t; A=sym([-21,19,-20; 19,-21,20; 40,-40,-40]);" + CR + "y0=[1; 0; -1]; y=expm(A*t)*y0")
End Sub
Sub g38()
CR = Chr(13) + Chr(10)
PutOnClipboard ("opt=odeset; opt.RelTol=1e-6;" + CR + "tic,[t,y]=ode45('c7xstf1',[0,1],[1;0;-1],opt); toc")
End Sub
Sub g39()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x1=exp(-2*t); x2=exp(-40*t).*cos(40*t); x3=exp(-40*t).*sin(40*t);" + CR + "y1=[0.5*x1+0.5*x2+0.5*x3, 0.5*x1-0.5*x2-0.5*x3, -x2+x3];" + CR + "plot(t,y,t,y1,':')")
End Sub
Sub g40()
CR = Chr(13) + Chr(10)
PutOnClipboard ("tic,[t2,y2]=rk_4('c7xstf1',[0,1,0.01],[1;0;-1]); toc")
End Sub
Sub g41()
CR = Chr(13) + Chr(10)
PutOnClipboard ("plot(t,y1,t2,y2,':')")
End Sub
Sub g42()
PutOnClipboard ("tic,[t2,y2]=ode45('c7exstf2',[0,100],[0;1]); toc")
End Sub
Sub g43()
PutOnClipboard ("length(t2), plot(t2,y2)")
End Sub
Sub g44()
CR = Chr(13) + Chr(10)
PutOnClipboard ("format long, [min(diff(t2)), max(diff(t2))]" + CR + "plot(t2(1:end-1),diff(t2))")
End Sub
Sub g45()
CR = Chr(13) + Chr(10)
PutOnClipboard ("opt=odeset; opt.RelTol=1e-6;" + CR + "tic,[t1,y1]=ode15s('c7exstf2',[0,100],[0;1],opt); toc")
End Sub
Sub g46()
PutOnClipboard ("length(t1), plot(t1,y1)")
End Sub
Sub g47()
CR = Chr(13) + Chr(10)
PutOnClipboard ("opt=odeset; opt.RelTol=1e-6;" + CR + "[t,x]=ode45('c7ximp',[0,10],[0; 0],opt); plot(t,x)")
End Sub
Sub g48()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[t,x]=ode15s('c7impode',[0,2],[1,0,0,1]); " + CR + "plot(t,x)")
End Sub
Sub g49()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=[1;0;0;1]; xd0=[0;0;1;1]; % 定义x0 和_x0" + CR + "res=ode15i('c7mimp',[0,2],x0,xd0); % 求解隐式微分方程" + CR + "plot(res.x,res.y) % 绘制各个状态的时间响应曲线")
End Sub
Sub g50()
CR = Chr(13) + Chr(10)
PutOnClipboard ("M=[1,0,0; 0,1,0; 0,0,0]; options=odeset; options.Mass=M;" + CR + "x0=[0.8; 0.1; 0.1]; " + CR + "[t,x]=ode15s('c7eqdae',[0,20],x0,options);" + CR + " plot(t,x)")
End Sub
Sub g51()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=[0.8; 0.1];" + CR + "fDae=inline(['[-0.2*x(1)+x(2)*(1-x(1)-x(2))+0.3*x(1)*x(2);','2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2))-2*x(2)*x(2)]'],'t','x');" + CR + "[t1,x1]=ode45(fDae,[0,20],x0); plot(t1,x1,t1,1-sum(x1'))")
End Sub
Sub g52()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x0=[0.8;0.1;2]; x0F=[1;1;0]; " + CR + "xd0=[1;1;1]; xd0F=[];" + CR + "[x0,xd0]=decic('c7midae',0,x0,x0F,xd0,xd0F);" + CR + " [x0,xd0] % 相容初始条件")
End Sub
Sub g53()
CR = Chr(13) + Chr(10)
PutOnClipboard ("res=ode15i('c7midae',[0,20],x0,xd0); " + CR + " plot(res.x,res.y)")
End Sub
Sub g54()
CR = Chr(13) + Chr(10)
PutOnClipboard ("options=odeset; options.Mass=@c7eqdae2; options.RelTol=1e-6;" + CR + "[t,x]=ode45(@c7fimp2,[0,10],[0;0],options); plot(t,x)")
End Sub
Sub g55()
CR = Chr(13) + Chr(10)
PutOnClipboard ("lags=[1 0.5];" + CR + " tx=dde23('c7exdde',lags,zeros(3,1),[0,10]);" + CR + "plot(tx.x,tx.y(2,:))")
End Sub
Sub g56()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[t,y]=shooting('c7fun1','c7fun2',[0,1],[1,2]);" + CR + " plot(t,y)")
End Sub
Sub g57()
CR = Chr(13) + Chr(10)
PutOnClipboard ("y0=((exp(2)-3)*exp(t)+(3-exp(1))*exp(2*t))/(4*exp(1)*(exp(1)-1))+3/4+t/2;" + CR + "norm(y(:,1)-y0)")
End Sub
Sub g58()
PutOnClipboard ("norm(y(end,1)-2)")
End Sub
Sub g59()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[t,y]=nlbound('c7fun4','c7fun3',[0,pi/2],[-1,1],1e-8);" + CR + " plot(t,y); " + CR + " set(gca,'xlim',[0,pi/2]);")
End Sub
Sub g60()
PutOnClipboard ("y0=tan(t-pi/4); norm(y(:,1)-y0)")
End Sub
Sub g61()
PutOnClipboard ("norm(y(end,1)-1)")
End Sub
Sub g62()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[t,y]=fdiff('c7fun5',[0,1],[1,4],50);" + CR + "  plot(t,y)")
End Sub
Sub g63()
CR = Chr(13) + Chr(10)
PutOnClipboard ("x=0:.05:1; t=0:0.05:2; m=0;" + CR + "sol=pdepe(m,@c7mpde,@c7mpic,@c7mpbc,x,t);" + CR + " surf(x,t,sol(:,:,1))")
End Sub
Sub g64()
CR = Chr(13) + Chr(10)
PutOnClipboard ("beta=8/3; rho=10; sigma=28;" + CR + "[t,x]=sim('c7mlor1b',[0,100]); plot(t,x)" + CR + "figure; plot3(x(:,1),x(:,2),x(:,3))")
End Sub
Sub g65()
CR = Chr(13) + Chr(10)
PutOnClipboard ("A1=[-13,3,-3; 106,-116,62; 207,-207,113];" + CR + "A2=diag([0.02,0.03,0.04]); B=[0; 1; 2];")
End Sub
Sub g66()
CR = Chr(13) + Chr(10)
PutOnClipboard ("[t,x]=sim('c7mdde3',[0,8]); plot(t,x)")
End Sub