基于MATLAB的M文件仿真

M文件：
k=1;
Int_F=inline('t','t');
for x=[1,3,5]
f_x(k)=x^3+x+log(x)*sin(x)+quad8(Int_F,0,x);
k=k+1;
end
f_x
>> Calcfx
Warning: QUAD8 is obsolete. We use QUADL instead. > In quad8 at 35
In Calcfx at 4
f_x =
2.5000 34.6550 140.9567
M文件：
function[mean,stdev]=stat(x)
n=length(x);
mean=sum(x)/n;
stdev=sqrt(sum(x-mean).^2/n);
>> x=[1,3,2];
>> [k,l]=stat(x)
k =
2
l =
微积分方程组的MA TLAB函数：
文件funcforex123.m
function xdot=funcforex123(t,x,flag,r,l,c)
xdot=zeros(2,1);
xdot(1)=-r/l*x(1)-1/l*x(2)+1/l*f(t);
xdot(2)=1/c*x(1);
function in=f(t)
in=(t>0)*1;
文件Ex123.m
l=1;
c=0.1;
for r=[1.5 3 5]
[t,x]=ode45('funcforex123',[-1,10],[0;0],[],r,l,c);
figure(1);plot(t,x(:,1));hold on;xlabel('time sec');
text(0.9,0.17,'\lefttarrow i_L(t)');grid;
figure(2);plot(t,x(:,2));hold on;xlabel('time sec');
text(0.5,0.3,'\leftarrow u_C(t)');grid;
End
>> ex123
Warning: Unable to interpret TeX string "\lefttarrow i_L(t)". > In ex123 at 5
Warning: Unable to interpret TeX string "\lefttarrow i_L(t)". > In ex123 at 7
Warning: Unable to interpret TeX string "\lefttarrow i_L(t)". > In ex123 at 7
文件ex123b.m
[t,x]=ode45('funcforex123',[-1,10],[0;0],[],2,1,0.1);
ts=0.001;
t1=-1:ts:10;
x1=interp1(t,x(:,2),t1,'spline');
plot(t1,x1,'k-.');
hold on;
x1dot=[diff([x1])/ts,0];
plot(t1,[diff([x1])/ts,0],'k');xlabel('time sec');
ht=10/sqrt(7.75).*exp(-1.5*t1).*sin(sqrt(7.75)*t1).*(t1>0);
plot(t1(1:50:length(t1)),ht(1:50:length(t1)),'ko');
legend('Syep response','Impulse response','Theoretic impluse response'); >> ex123b
欧拉算法的MA TLAB程序：
文件sybeuler.m
function [tout,yout]=sybeuler(odefile,t0,h,th,y0,P) tout=[t0:h:th]';
yout(length(tout),length(y0))=0;
kk=1;
for t=tout'
yout(kk,:)=y0';
kk=kk+1;
k1=h*eval([odefile'(t,y0,P)']);
y0=y0+k1;
end
文件mystateEQ.m
function xdot=mystateEQ(t,x,P)
xdot=zeros(2,1);
r=P(1);l=P(2);c=P(3);
xdot(1)=-r/l*x(1)-1/l*x(2)+1/l*f(t);
xdot(2)=1/c*x(1);
function in=f(t)
in=(t>0)*1;
文件mystateEQforODE45.m
function xdot=mystateEQforODE45(t,x,flag,P) xdot=zeros(2,1);
r=P(1);l=P(2);c=P(3);
xdot(1)=-r/l*x(1)-1/l*x(2)+1/l*f(t);
xdot(2)=1/c*x(1);
function in=f(t)
in=(t>0)*1;
文件sybeulerTEST.m
tic;
[t,x]=ode45('mystateEQforODE45',[0,10],[0;0],[],[3 1 0.1]); toc
plot(t,x,'k.');hold on;
for h=[0.2,0.1,0.01]
tic;
[t,x]=sybeuler('mystateEQ',0,h,10,[0;0],[3 1 0.1]);
toc
plot(t,x,'k');
end
legend('u_C(t)ode45','i_L(t)ode45','sybeuler');。