重要的程序

用matlab解延迟微分方程(四介Rungkuter)： 下面是个例子： % 欲解方程组： % dx/dt=y-x+x(t-1); % dy/dt=y-xy+y(t-1); x0=-1; %初试条件 y0=1; %初试条件 t=0:01:3; yy=zeros(size(t)); xx=zeros(size(t)); yy(1)=x0; yy(1)=y0; dt=t(2)-t(1); Nk=round(1/dt); for k=1:length(t)-1;     if t(k)<1;         k1=yy(k)-xx(k);         g1=yy(k)-xx(k)*yy(k);                k2=yy(k)+0.5*g1*dt-[xx(k)+0.5*dt*k1];         g2=yy(k)+0.5*g1*dt-[xx(k)+0.5*dt*k1]*[yy(k)+0.5*g1*dt];         k3=yy(k)+0.5*g2*dt-[xx(k)+0.5*dt*k2];         g3=yy(k)+0.5*g2*dt-[xx(k)+0.5*dt*k1]*[yy(k)+0.5*g2*dt];         k4=yy(k)+g3*dt-[xx(k)+dt*k3];         g4=yy(k)+g3*dt-[xx(k)+dt*k3]*[yy(k)+g3*dt];     else         k1=yy(k)-xx(k)+xx(k-Nk);         g1=yy(k)-xx(k)*yy(k)+yy(k-Nk);                k2=yy(k)+0.5*g1*dt-[xx(k)+0.5*dt*k1]+xx(k-Nk)+0.5*dt*k1;         g2=yy(k)+0.5*g1*dt-[xx(k)+0.5*dt*k1]*[yy(k)+0.5*g1*dt]+yy(k-Nk)+0.5*dt*g1;         k3=yy(k)+0.5*g2*dt-[xx(k)+0.5*dt*k2]+xx(k-Nk)+0.5*dt*k2;         g3=yy(k)+0.5*g2*dt-[xx(k)+0.5*dt*k1]*[yy(k)+0.5*g2*dt]+yy(k-Nk)+0.5*dt*g2;         k4=yy(k)+g3*dt-[xx(k)+dt*k3]+xx(k-Nk)+dt*k3;         g4=yy(k)+g3*dt-[xx(k)+dt*k3]*[yy(k)+g3*dt]+yy(k-Nk)+dt*k3;            end     xx(k+1)=xx(k)+dt*[k1+2*k2+2*k3+k4]/6;     yy(k+1)=yy(k)+dt*[g1+2*g2+2*g3+g4]/6; end subplot(121);plot(t,xx);title('t~x(t)'); subplot(122);plot(t,yy);title('t~y(t)'); 4阶变步长龙格-库塔法解dy=fun(x,y)型方程 function change_step_RK(fun); % 4阶变步长龙格-库塔法解dy=fun(x,y)型方程 % Example: %     fun=inline('cos(x)+sin(y)'); %     change_step_RK(fun); AbsTol=1e-4; h=0.01; p=4;   % p是对应的阶数 x0=0;  % x0是起点 xN=1;  % xN是终点 y0=0;  % 初值 x=x0; y=y0; n=1; p21=2^p-1; while x(end)AbsTol;            x1=x2;            y1=y2;            h=h/2;            [x2,y2]=RK_f(fun,x(n),y(n),h/2);        end    end    [xa,ya]=RK_f(fun,h,x(n),y(n));    x(n+1)=xa;    y(n+1)=ya;    n=n+1; end plot(x,y,'k'); function [xa,ya]=RK_f(fun,h,x,y); % 4阶龙格-库塔法解下一点的值 k1=fun(x,y); k2=fun(x+h/2,y+h*k1/2); k3=fun(x+h/2,y+h*k2/2); k4=fun(x+h,y+h*k3); xa=x+h; ya=y+h*(k1+k2*2+2*k3+k4)/6; 脉冲作用下的解微分方程matlab求解一例： Differential equations: dx1/dt=-x1/t1+x2+x3*(x2-x1) dx2/dt=r(t)-x2-x3*(x2-x1) dx3/dt=A*x2+x3*(x2-x1)-x3/t0 initial conditions: x1(0)=x2(0)=x3(0)=0; r(t) is as follow: mr=mod(t,tp+ti); r=0; if mr<2;    r=4; end here is : tp=2;ti=50;t0=10;t1=0.5;A=1.5e-8; range of time: from 0 to 200. the program is as follow: function dx=cyj(t,x); % 脉冲作用下的解微分方程 % Example: %    [t,y]=ode45('cyj',[0,200],[0,0,0]); %    plot(t,y(:,1),t,y(:,2),t,y(:,3)); %    legend('x1','x2','x3',0); % tp=2; ti=50; mr=mod(t,tp+ti); r=0; if mr<2;    r=4; end t0=10; t1=0.5; A=1.5e-8; Ji=x(3)*(x(2)-x(1)); dx=[-x(1)/t1+x(2)+Ji;r-x(2)-Ji;A*x(2)+Ji-x(3)/t0]; matlab解多元微分方程组一例： 问题[by 半月]： 具体题目如下: 令mj>0 (j=0,1,2)为三个物体的质量,xj∈R^2和xj'∈R^2是在惯性坐标系下它们的位置和速度(列)向量.用rj

表

关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf

示两个物体之间的距离 r0=|x2-x1|, r1=|x0-x2|, r2=|x1-x0|, 如果选择长度,时间和质量单位,使得引力常数等于1,则牛顿运动方程为 x0''=m1*(x1-x0)/r2^3+m2*(x2-x0)/r1^3 x1''=m2*(x2-x1)/r0^3+m0*(x0-x1)/r2^3 x2''=m0*(x0-x2)/r0^3+m1*(x1-x2)/r0^3 对给定的初始条件: m0=5, m1=3, m2=4; x0=[1;-1], x1=[1;3], x2=[-2;-1]; x0'=x1'=x2'=0 画出质点运动的轨迹,时间为t(0x0 r1=abs(y(1)-y(3));   % y(2)--->x1 r2=abs(y(2)-y(1));   % y(3)--->x2 dy(1,1)=y(4); dy(2,1)=y(5); dy(3,1)=y(6); dy(4,1)=m1*(y(2)-y(1))/r2^3+m2*(y(3)-y(1))/r1^3;   dy(5,1)=m2*(y(3)-y(2))/r0^3+m0*(y(1)-y(2))/r2^3; dy(6,1)=m0*(y(1)-y(3))/r0^3+m1*(y(2)-y(3))/r0^3; % x0''=m1*(x1-x0)/r2^3+m2*(x2-x0)/r1^3 % x1''=m2*(x2-x1)/r0^3+m0*(x0-x1)/r2^3 % x2''=m0*(x0-x2)/r0^3+m1*(x1-x2)/r0^3 又一例程[for zyang]： function dx=zyang(t,x); % Example: %   x0=[0.3,1.5,0.2,0.95]; %   [t,x]=ode45('zyang',[0,100],x0); %   plot(x(:,2),x(:,4)); %   figure;plot(t,x(:,4)); q1=10.82; dx(1)=x(2); dx(2)=16.6667*sin(x(3)-x(1)+0.0873)*x(4)-...    0.1667*x(2)+1.8807; dx(3)=496.8718*x(4)*x(4)-166.6667*cos(x(3)-...    x(1)-0.0873)*x(4)-666.6667*cos(x(3)-...    0.2094)*x(4)-93.3333*x(4)+33.3333*q1+43.3333; dx(4)=-78.7638*x(4)*x(4)+26.2172*cos(x(3)-...    x(1)-0.0124)*x(4)+104.8689*cos(x(3)-...    0.1346)*x(4)+14.5229*x(4)-5.2288*q1-7.0327; dx=dx'; 不同inline定义的函数之间的调用： 下面的例子，实现Df对gx的调用 B=1; m0=-1/7; m1=-40/7; clc; gx=inline('m0*x+0.5*(m1-m0)*abs(abs(x+B)-abs(x-B))',...    'x','B','m0','m1'); a=gx(1,B,m0,m1) Df=inline('2*gx(x,B,m0,m1)','x','gx','B','m0','m1') 关于求解变参数lorenz方程的程序： 怎么求解变参数lorenz方程，比如，改变参数r，使r的取值在0.1到30之间，步长去为0.1，将结果放入一个细胞数组中。直接用ode45好像不能向里面传递参数，还有什么其他的方法解决么？ 程序如下: % 关于求解变参数lorenz方程 fun=inline(['[-b*x(1)+x(2)*x(3);',...            '-xigma*(x(2)-x(3));',...            'r*x(2)-x(3)-x(1)*x(2)]'],...            't','x','flag','xigma','r','b'); xigma=10;b=8/3; rr=0.1:.1:3; Data_T=cell(length(rr),1); Data_X=cell(length(rr),1); tspan=[0,10]; x0=[0,0,1e-10]; for k=1:length(rr);    r=rr(k);    [T,X]=ode45(fun,tspan,[0,0,1],[],xigma,r,b);    Data_T(k)={T};    Data_X(k)={X}; end p=randperm(length(rr));p(2:end)=[]; plot(cell2mat(Data_T(1)),cell2mat(Data_X(1))); title(['p=',num2str(p)]); A=cell2mat(Data_X(p)); figure;plot3(A(:,1),A(:,2),A(:,3)); title(['p=',num2str(p)]); 一个微分方程解的动态显示： fun=inline('[y(2);a*y(1)-2*x*y(2)]','x','y','flag','a'); tspan=[0,20]; y0=[0,1]; g1=plot(0,0.5);hold on; g2=plot(0,0.6); xlim(tspan) for a=-2:.1:2;    [xq,yq]=ode45(fun,tspan,y0,[],a);    set(g1,'xdata',xq,'ydata',yq(:,1));    set(g2,'xdata',xq,'ydata',yq(:,2));    title(['\lambda = ',num2str(a)]);    pause(0.5); end 基于matlab的电路仿真：close all;clear;clc; figure('position',[189 89 714 485]); Na=['输入波形[请选择]|输入波形:正弦波|',...         '输入波形:方形波|输入波形:脉冲波']; Ns={'sin','square','pulse'};R=2;  % default parameters: resistance C=2;  % default parameters: capacitance f=10; % default parameters: frequencyTAU=R*C;; tff=10;  % length of time ts=1/f;  % sampling lengthsys1=tf([1],[1,1]);     % systems for integral circuit sys2=tf([1,0],[1,1]);   % systems for differential circuita1=axes('position',[0.1,0.6,0.3,0.3]); po1=uicontrol(gcf,'style','popupmenu',...     'unit','normalized','position',[0.15,0.9,0.2,0.08],...     'string',Na,'fontsize',12,'callback',[]); set(po1,'callback',['KK=get(po1,''Value'');if KK>1;',...         'st=char(Ns(KK-1));[U,T]=gensig(st,R*C,tff,1/f);',...         'axes(a1);plot(T,U);ylim([min(U)-0.5,max(U)+0.5]);',...         'end;']); Ma=['电路类型[请选择]|电路类型:积分型|电路类型:微分型']; a2=axes('position',[0.5,0.6,0.3,0.3]);box on; set(gca,'xtick',[]);set(gca,'ytick',[]); po2=uicontrol(gcf,'style','popupmenu',...     'unit','normalized','position',[0.55,0.9,0.2,0.08],...     'string',Ma,'fontsize',12,'callback',[]); set(po2,'callback',['KQ=get(po2,''Value'');axes(a2);',...         'if KQ==1;cla;elseif KQ==2;',...         'plot(0.14+0.8i+0.02*exp(i*[0:.02:8]),''k'');hold on;',...         'plot(0.14+0.2i+0.02*exp(i*[0:.02:8]),''k'');',...         'plot(0.84+0.2i+0.02*exp(i*[0:.02:8]),''k'');',...         'plot(0.84+0.8i+0.02*exp(i*[0:.02:8]),''k'');',...         'plot([0.16,0.82],[0.2,0.2],''k'');',...         'plot([0.16,0.3],[0.8,0.8],''k'');',...         'plot([3,4,4,3,3]/10,[76,76,84,84,76]/100,''k'');',...         'plot([0.4,0.82],[0.8,0.8],''k'');',...         'plot([0.6,0.6],[0.8,0.53],''k'');',...         'plot([0.6,0.6],[0.2,0.48],''k'');',...         'plot([0.55,0.65],[0.53,0.53],''k'');',...         'plot([0.55,0.65],[0.48,0.48],''k'');',...         'text(0.33,0.7,''R'');',...         'text(0.63,0.42,''C'');',...         'text(0.8,0.5,''U_o'');text(0.1,0.5,''U_i'');hold off;',...         'set(gca,''xtick'',[]);set(gca,''ytick'',[]);',...         'else;',...         'plot(0.14+0.8i+0.02*exp(i*[0:.02:8]),''k'');hold on;',...         'plot(0.14+0.2i+0.02*exp(i*[0:.02:8]),''k'');',...         'plot(0.84+0.2i+0.02*exp(i*[0:.02:8]),''k'');',...         'plot(0.84+0.8i+0.02*exp(i*[0:.02:8]),''k'');',...         'plot([0.16,0.82],[0.2,0.2],''k'');',...         'plot([0.16,0.33],[0.8,0.8],''k'');',...         'plot([0.33,0.33],[0.74,0.86],''k'');',...         'plot([0.36,0.36],[0.74,0.86],''k'');',...         'plot([0.36,0.82],[0.8,0.8],''k'');',...         'plot([0.6,0.6],[0.8,0.58],''k'');',...         'plot([0.6,0.6],[0.2,0.42],''k'');',...         'plot([58,62,62,58,58]/100,[42,42,58,58,42]/100,''k'');',...         'text(0.33,0.68,''C'');',...         'text(0.63,0.42,''R'');',...         'text(0.8,0.5,''U_o'');text(0.1,0.5,''U_i'');hold off;',...         'set(gca,''xtick'',[]);set(gca,''ytick'',[]);',...         'end;axis([0,1,0,1]);']); a3=axes('position',[0.1,0.1,0.3,0.3]); pu3=uicontrol(gcf,'style','push',...     'unit','normalized','position',[0.15,0.43,0.2,0.04],...     'string','显示输出波形','fontsize',12,...     'callback',['axes(a3);if KQ==2;lsim(sys1,U,T);',...         'elseif KQ==3;lsim(sys2,U,T);else,cla;end;']); uicontrol(gcf,'style','frame',...     'unit','normalized','position',[0.47,0.13,0.5,0.32],...     'string','10','BackgroundColor',[0.8,0.8,0.8],...     'fontsize',16,'fontname','times new roman'); s1=uicontrol(gcf,'style','slider',...     'unit','normalized','position',[0.55,0.38,0.3,0.03],...     'Value',0.5,'SliderStep',[0.01 0.01],...     'callback',['R=1+get(s1,''Value'')*2;',...         '[U,T]=gensig(st,R*C,tff,1/f);',...         'set(t1,''string'',num2str(R));']); uicontrol(gcf,'style','text',...     'unit','normalized','position',[0.5,0.37,0.05,0.05],...     'string','R','BackgroundColor',[0.8,0.8,0.8],...     'fontsize',16,'fontname','times new roman'); t1=uicontrol(gcf,'style','text',...     'unit','normalized','position',[0.86,0.37,0.1,0.05],...     'string','2','BackgroundColor',[0.8,0.8,0.8],...     'fontsize',16,'fontname','times new roman'); s2=uicontrol(gcf,'style','slider',...     'unit','normalized','position',[0.55,0.28,0.3,0.03],...     'Value',0.5,'SliderStep',[0.01 0.01],...     'callback',['C=1+get(s2,''Value'')*2;',...         '[U,T]=gensig(st,R*C,tff,1/f);',...         'set(t2,''string'',num2str(C));']); uicontrol(gcf,'style','text',...     'unit','normalized','position',[0.5,0.27,0.05,0.05],...     'string','C','BackgroundColor',[0.8,0.8,0.8],...     'fontsize',16,'fontname','times new roman'); t2=uicontrol(gcf,'style','text',...     'unit','normalized','position',[0.86,0.27,0.1,0.05],...     'string','2','BackgroundColor',[0.8,0.8,0.8],...     'fontsize',16,'fontname','times new roman'); set(gcf,'NumberTitle','off','MenuBar','none'); set(gcf,'name','基于matlab的电路仿真'); uu=uimenu(gcf,'label','About'); uimenu(uu,'label','Author:Zjliu'); ku=uimenu(gcf,'label','Readme'); kw1=uimenu(gcf,'label','执行顺序:'); kw2=uimenu(gcf,'label','设定参数-->确定波形-->选择电路类型-->输出波形'); set(kw1,'ForegroundColor','r');set(kw2,'ForegroundColor','b'); uicontrol(gcf,'style','text',...     'unit','normalized','position',[0.5,0.47,0.46,0.05],...     'string','请选择波形之前改变下面的参数',...     'BackgroundColor',[0.8,0.8,0.8],...     'fontsize',16,'fontname','times new roman');    % readme s3=uicontrol(gcf,'style','slider',...     'unit','normalized','position',[0.55,0.18,0.3,0.03],...     'Value',0.5,'SliderStep',[0.01 0.01],...     'callback',['f=5+get(s3,''Value'')*10;',...         '[U,T]=gensig(st,R*C,tff,1/f);',...         'set(t3,''string'',num2str(f));']); uicontrol(gcf,'style','text',...     'unit','normalized','position',[0.5,0.17,0.05,0.05],...     'string','f','BackgroundColor',[0.8,0.8,0.8],...     'fontsize',16,'fontname','times new roman'); t3=uicontrol(gcf,'style','text',...     'unit','normalized','position',[0.86,0.17,0.1,0.05],...     'string','10','BackgroundColor',[0.8,0.8,0.8],...     'fontsize',16,'fontname','times new roman');