matlab ode45和矩阵生成有向网络图

Matlab中解常微分方程的ode45
ode是专门用于解微分方程的功能函数，他有ode23,ode45,ode23s等等，采用的是Runge-Kutta算法。

ode45表示采用四阶，五阶runge-kutta单步算法,截断误差为(Δx)^3。

解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,换用ode23来解.其他几个也是类似的用法
使用方法
[T,Y] = ode45(odefun,tspan,y0)
odefun 是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名
tspan 是区间
[t0 tf] 或者一系列散点[t0,t1,...,tf]
y0 是初始值向量
T 返回列向量的时间点
Y 返回对应T的求解列向量
[T,Y] = ode45(odefun,tspan,y0,options)
options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等
[T,Y,TE,YE,IE] =ode45(odefun,tspan,y0,options)
每组(t,Y)之产生称为事件函数。

每次均会检查是否函数等于零。

并决定是否在零时终止运算。

这可以在函数中之特性上设定。

例如以events 或@events产生一函数。

[value, isterminal,direction]=events(t,y)其中，value(i)为函数之值，isterminal(i)=1时运算在等于零时停止，=0时继续；direction(i)=0时所有零时均需计算（默认值），+1在事件函数增加时等于零，-1在事件函数减少时等于零等状况。

此外，TE, YE, IE则分别为事件发生之时间，事件发生时之答案及事件函数消失时之指针i。

sol =ode45(odefun,[t0 tf],y0...) sol 结构体输出结果
应用举例
1 求解一阶常微分方程
程序:
) (y+3*t)/t^2; %定义函数
tspan=[1 4]; %求解区间
y0=-2; %初值
[t,y]=ode45(odefun,tspan,y0);
plot(t,y) %作图
title('t^2y''=y+3t,y(1)=-2,1<t<4')
legend('t^2y''=y+3t') xlabel('t')
ylabel('y') % 精确解
% dsolve('t^2*Dy=y+3*t','y(1)=-2')
% ans =
% (3*Ei(1) - 2*exp(1))/exp(1/t) - (3*Ei(1/t))/exp(1/t)
2 求解高阶常微分方程
关键是将高阶转为一阶,odefun的书写.
F(y,y',y''...y(n-1),t)=0用变量替换,y1=y,y2=y'...注意odefun方程定义为列向量dxdy=[y(1),y(2)....]
程序:
function Testode45
tspan=[3.9 4.0]; %求解区间
y0=[2 8]; %初值
[t,x]=ode45(@odefun,tspan,y0);
plot(t,x(:,1),'-o',t,x(:,2),'-*')
legend('y1','y2')
title('y'' ''=-t*y + e^t*y'' +3sin2t')
xlabel('t') ylabel('y')
function y=odefun(t,x)
y=zeros(2,1); % 列向量
y(1)=x(2);
y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);
end
end
matlab练习程序（矩阵生成有向网络图）
早知道有向图和无向图差别没有想象中的大我就写到一起了。

函数中使用的arrow画箭头函数是在这个网站下的。

%函数名netplot
%使用方法输入请help netplot
%无返回值
%函数只能处理有向图
%作者：tiandsp
%最后修改：2012.12.26
function netplot(A,flag)
%调用方法输入netplot(A,flag)，无返回值
%A为邻接矩阵或关联矩阵
%flag=1时处理邻接矩阵
%flag=2时处理关联矩阵
%函数只能处理有向图
if flag==1 %邻接矩阵表示有向图
D_netplot(A);
return;
end
if flag==2 %关联矩阵表示有向图
[m n]=size(A); %关联矩阵变邻接矩阵
W=zeros(m,m);
for i=1:n
a=find(A(:,i)~=0);
if A(a(1),i)==1 W(a(1),a(2))=1;
else W(a(2),a(1))=1;
end
end D_netplot(W);
return;
end
function D_netplot(A)
[n n]=size(A);
w=floor(sqrt(n));
h=floor(n/w);
x=[];
y=[];
for i=1:h %使产生的随机点有其范围，使显示分布的更广
for j=1:w
x=[x 10*rand(1)+(j-1)*10];
y=[y 10*rand(1)+(i-1)*10];
end
end ed=n-h*w;
for i=1:ed
x=[x 10*rand(1)+(i-1)*10];
y=[y 10*rand(1)+h*10];
end plot(x,y,'r*');
title('网络拓扑图');
for i=1:n
for j=1:n
if A(i,j)~=0
c=num2str(A(i,j)); %将A中的权值转化为字符型
text((x(i)+x(j))/2,(y(i)+y(j))/2,c,'Fontsize',10); %显示边的权值
arrow([x(i) y(i)],[x(j) y(j)]); %带箭头的连线
end
text(x(i),y(i),num2str(i),'Fontsize',14,'color','r'); %显示点的序号 hold on;
end
end
end
end
运行结果：。