有關微分方程編輯 runge-kutta 的問題

eatnche594

這是我的副程式

function  yp=ivp13(t,y)
global a b c d
yp=zeros(1,14);  % note: column vector
yp(1)=y(2);
yp(2)=((6.15*10.^6*(y(3)-y(1)))-(6.15*10.^6*y(1)))/(17.88336)-yp(14);
yp(3)=y(4);
yp(4)=((8.1*10.^7*(y(5)-y(3)))-(6.15*10.^6*(y(3)-y(1))))/(0.088062)-yp(14);
yp(5)=y(6);
yp(6)=(6.03*10.^6*(y(7)-y(5))-8.1*10.^7*(y(5)-y(3)))/(0.1544472)-yp(14);
yp(7)=y(8);
yp(8)=((b+0.5*c*y(7)+0.5*d*y(7).^2)*2.01*10.^6*(a+y(9)-(a+b*y(7)+0.5*c*y(7).^2)))/(0.0758688)-yp(14);
yp(9)=y(10);
yp(10)=(2.01*10.^6*(y(11)-y(9))-2.83*10.^8.*(a+y(9)-(a+b*y(7)+0.5*c*y(7).^2)))/(0.023709)-yp(14)*b-c*y(14).^2;
yp(11)=y(12);
yp(12)=(1000-2.01*10.^6*(y(11)-y(9)))/(0.3387)-yp(14)*b-c*y(14).^2;
yp(13)=y(14);
yp(14)=(107.257-a*(0.25+0.104+5)*(b*y(14).^2+1000))/(0.096+13.2+1.3+2.28+1.12+((0.25+0.104+5)*a.^2));
yp=yp';

這執行檔
clc;clear all;
load data_a.txt; load data_b.txt; load data_c.txt; load data_d.txt
tspan=[0 0.5]; x0=[0;0;0;0;0;0;0;0;0;0;0;0;0;0];
global a b c d
for k=1:58
   a=data_a(k); b=data_b(k); c=data_c(k); d=data_d(k);
   [ta,xa]=ode45('ivp1',tspan,x0);
end
for k=59:303
   a1=data_a(k); b1=data_b(k); c1=data_c(k); d1=data_d(k);
   [ta,xa]=ode45('ivp1',tspan,x0);
end
for k=303:359
   a2=data_a(k); b2=data_b(k); c2=data_c(k); d2=data_d(k);
   [ta,xa]=ode45('ivp1',tspan,x0);
end
plot(ta,xa(:,12))

當我將tspan改的越大,為何跑程式的時間會越長而且是跑3~6天
超過0.6 時 就會出現 out of memory
有那位大大 可以幫幫我看看 那裏可以改或者 哪裡錯了?
附件為data_a.txt 、data_b.txt、 data_c.txt、 data_d.txt

[ 本帖最后由 ChaChing 于 2009-4-17 08:35 编辑 ]