| 网站首页 | 业界新闻 | 小组 | 威客 | 人才 | 下载频道 | 博客 | 代码贴 | 在线编程 | 编程论坛
欢迎加入我们,一同切磋技术
用户名:   
 
密 码:  
共有 5595 人关注过本帖
标题:求助,matlab求解微分方程组
只看楼主 加入收藏
mon890
Rank: 1
等 级:新手上路
帖 子:4
专家分:0
注 册:2006-9-15
收藏
 问题点数:0 回复次数:6 
求助,matlab求解微分方程组

求助,很着急,求matlab高手帮帮忙吧

我想解一个微分方程组,主程序是这样的
a=0.000000012;
b=0.000000024;
M=3000;
h=(b-a)/M;
t=zeros(1,M+1);
x=zeros(1,M+1);
y=zeros(1,M+1);
z=zeros(1,M+1);
t=a:h:b;
x(1)=6.4164e+010;
y(1)=4372.7;
z(1)=2.2129e+024;
for j=1:M

kx1=h*feval('fx',t(j),x(j),y(j),z(j));
ky1=h*feval('fy',t(j),x(j),y(j),z(j));
kz1=h*feval('fz',t(j),x(j),y(j));


kx2=h*feval('fx',t(j)+h/2,x(j)+kx1/2,y(j)+ky1/2,z(j)+kz1/2);
ky2=h*feval('fy',t(j)+h/2,x(j)+kx1/2,y(j)+ky1/2,z(j)+kz1/2);
kz2=h*feval('fz',t(j)+h/2,x(j)+kx1/2,z(j)+kz1/2);


kx3=h*feval('fx',t(j)+h/2,x(j)+kx2/2,y(j)+ky2/2,z(j)+kz2/2);
ky3=h*feval('fy',t(j)+h/2,x(j)+kx2/2,y(j)+ky2/2,z(j)+kz2/2);
kz3=h*feval('fz',t(j)+h/2,x(j)+kx2/2,z(j)+kz2/2);


kx4=h*feval('fx',t(j)+h,x(j)+kx3,y(j)+ky3,z(j)+kz3);
ky4=h*feval('fy',t(j)+h,x(j)+kx3,y(j)+ky3,z(j)+kz3);
kz4=h*feval('fz',t(j)+h,x(j)+kx3,z(j)+kz3);


x(j+1)=x(j)+(kx1+2*kx2+2*kx3+kx4)/6;
y(j+1)=y(j)+(ky1+2*ky2+2*ky3+ky4)/6;
z(j+1)=z(j)+(kz1+2*kz2+2*kz3+kz4)/6;
end

plot(t,x);

调用的fx,fy,fz分别为
function dx=fx(t,x,y,z)
Gn=8.4e-13;nth=2.0e+24;k=0.015;tin=8e-12;
dx=(1/2)*Gn*(z-nth)*x+(k/tin)*C1*cos(y-C2);
function dy=fy(t,x,y,z)
Gn=8.4e-13;n0=1.4e+24;tin=8e-12;k=0.015;a=6;
dy=(1/2)*a*Gn*(z-n0)-(k/tin)*(C1/x)*sin(y-C2);
function dz=fz(t,x,z)
Gn=8.4e-13;e=1.6e-19;V=1.0e-16;ts=2e-9;Ith=54e-3;n0=1.4e+24;
dz=(1.3*Ith)/(e*V)-(1/ts)*z-Gn*(z-n0)*x^2;

上边两个式子中有C1,C2,是两个数组,分别为

C1=[8.2e-011, -1.4679e-007, 2.2247e-007, 2.3087e-007, -1.0047e-007, 3.153e-008, 3.2047e-007, 1.0044e-007, -1.7736e-007, 1.8828e-007, -1.9584e-007, -1.1957e-007, 2.6189e-007, 2.4234e-008, 1.6046e-007, -1.2331e-007, -1.9984e-007, 4.6585e-007, -2.9163e-007, 1.4481e-008, -1.2686e-007, 6.6135e-007, -1.1262e-007, 1.3977e-007, -6.5892e-008, -3.6625e-007, -2.467e-007, 7.4485e-008, 1.475e-007, 2.7209e-007, 2.4435e-007, 1.0191e-007, 8.4376e-008, 1.5446e-008, -1.4789e-007, -1.9714e-007, 3.0642e-007, -3.1694e-007, 1.3063e-007, 1.2397e-007, -5.0877e-008, 8.2291e-008, -4.9396e-008, 4.4198e-008, -6.6902e-010, -2.0217e-008, 1.4079e-008, 2.1972e-008, 7.8652e-009, 1.0089e-008, 2.0888e-008, -5.0502e-008, 1.75e-007, 8.8072e-007, 1.6731e-006, 2.1356e-006, 2.0616e-006, 1.4339e-006, 2.6416e-007, -1.1631e-006, -2.3767e-006, -2.8328e-006, -1.3912e-006, 3.3939e-006, 1.7953e-005, 0.00010396, 0.00069155, 0.0052414, 0.045308, 0.44673, 5.0174, 64.156, 933.44, 15451, 2.9058e+005, 6.2072e+006, 1.5046e+008, 4.1381e+009, 1.2671e+011, 4.4426e+011, 8.9674e+010, 1.8253e+010, 4.2933e+009, 1.1719e+009, 3.7112e+008, 1.3615e+008, 5.7867e+007, 2.8467e+007, 1.6195e+007, 1.0647e+007, 8.0882e+006, 7.0847e+006, 7.1595e+006, 8.3524e+006, 1.1181e+007, 1.7265e+007, 3.0646e+007, 6.2556e+007, 1.4634e+008, 3.9367e+008, 1.2152e+009, 4.2877e+009, 1.7323e+010, 7.8254e+010, 2.8099e+011, 2.251e+011, 8.7931e+010, 3.4687e+010, 1.5524e+010, 8.0087e+009, 4.7717e+009, 3.2836e+009, 2.609e+009, 2.391e+009, 2.5263e+009, 3.0734e+009, 4.3131e+009, 6.934e+009, 1.2876e+010, 2.7086e+010, 6.4164e+010,
];
C2=[0, -5.3281, -9.5561, -12.689, -14.733, -15.694, -15.576, -14.385, -12.126, -8.8054, -4.4275, 1.0022, 7.4783, 14.996, 23.549, 33.133, 43.744, 55.374, 68.021, 81.677, 96.34, 112, 128.66, 146.31, 164.94, 184.56, 205.15, 226.71, 249.24, 272.73, 297.17, 322.57, 348.91, 376.2, 404.42, 433.58, 463.66, 494.67, 526.59, 559.44, 593.19, 627.84, 663.4, 699.85, 737.2, 775.43, 814.54, 854.54, 895.4, 937.14, 979.74, 1023.2, 1067.5, 1112.7, 1158.7, 1205.6, 1253.3, 1301.8, 1351.2, 1401.4, 1452.4, 1504.3, 1556.9, 1610.4, 1664.7, 1719.7, 1775.6, 1832.3, 1889.8, 1948, 2007, 2066.8, 2127.4, 2188.8, 2250.9, 2313.8, 2377.5, 2441.9, 2506.9, 2558.9, 2593.8, 2628.8, 2664.6, 2701.4, 2739, 2777.5, 2816.9, 2857.1, 2898.3, 2940.3, 2983.2, 3026.9, 3071.5, 3116.9, 3163.2, 3210.3, 3258.3, 3307.1, 3356.7, 3407.1, 3458.4, 3510.5, 3563.4, 3617, 3669.2, 3712.4, 3751.3, 3790.2, 3829.9, 3870.5, 3911.9, 3954.2, 3997.3, 4041.3, 4086.2, 4131.9, 4178.4, 4225.8, 4274, 4323, 4372.7,
];
我的意思就是想在主程序里,那个for 循环,每循环一次,就引用数组里的一个数值,主程序里算出来的数再接着继续存储在这两个数组里以便能继续算下去。拜托哪位高手帮帮我,我已经请教了很多人了,可是没弄出来,希望能有人帮忙,如果还有不清楚地,可以回个帖子,或者加我QQ也行,我的QQ:27866661,很着急,希望斑竹帮帮忙,推荐一下


搜索更多相关主题的帖子: matlab 微分方程组 feval zeros 求解 
2006-11-06 16:17
hitzhang
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:版主
威 望:21
帖 子:369
专家分:52
注 册:2006-9-24
收藏
得分:0 
你的t向量是0阿!调用的fx,fy,fz里似乎没有t阿
另外C1,C2是怎么算出来的
3000次循环而他们的维数只有121,是没算完吗?

如果每循环一次,就引用数组里的一个数值,主程序里算出来的数再接着继续存储在这两个数组里以便能继续算下去
那么可以把fx,fy,fz添加一个输入参量C

2006-11-07 14:17
mon890
Rank: 1
等 级:新手上路
帖 子:4
专家分:0
注 册:2006-9-15
收藏
得分:0 
C1C2是另外的程序已经算出来的,怎么添加输入参量啊,这样说我也不太明白,能在我那个程序上直接改吗?谢谢啦!
2006-11-08 18:11
hitzhang
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:版主
威 望:21
帖 子:369
专家分:52
注 册:2006-9-24
收藏
得分:0 

调用的fx,fy,fz里似乎没有t阿
function dx=fx(t,x,y,z)
Gn=8.4e-13;nth=2.0e+24;k=0.015;tin=8e-12;
dx=(1/2)*Gn*(z-nth)*x+(k/tin)*C1*cos(y-C2);
c1,c2是数组算出来的dx也是数组,我不太明白你的具体问题
可以这样添加参量
function dx=fx(t,x,y,z,C1,C2)
Gn=8.4e-13;nth=2.0e+24;k=0.015;tin=8e-12;
dx=(1/2)*Gn*(z-nth)*x+(k/tin)*C1*cos(y-C2);
在主程序中
for j=1:M

kx1=h*feval('fx',t(j),x(j),y(j),z(j),C1(j),C2(j));
.
.
.


2006-11-11 11:15
mon890
Rank: 1
等 级:新手上路
帖 子:4
专家分:0
注 册:2006-9-15
收藏
得分:0 
我自己编的那个可能有点毛病,这样吧,你能帮我看看能不能编个程序,解这个微分方程组
dx(t)/dt=(1/2)*Gn*(z(t)-nth)*x(t)+(k/tin)*x(t-tt)*cos(y(t)-y(t-tt));
dy(t)/dt=(1/2)*a*Gn*(z(t)-n0)-(k/tin)*(x(t-tt)/x(t))*sin(y(t)-y(t-tt));
dz(t)/dt=(1.3*Ith)/(e*V)-(1/ts)*z(t)-Gn*(z(t)-n0)*x(t)^2;
Gn=8.4e-13;nth=2.0e+24;k=0.015;tin=8e-12;
n0=1.4e+24;k=0.015;a=6;
e=1.6e-19;V=1.0e-16;ts=2e-9;Ith=54e-3;tt=1.2e-9;
我只知道 用4阶龙格库塔法好像能解,别的我也不知道了,拜托斑竹一定帮我看看啊,已经找了好多人了,都没弄明白,郁闷阿
2006-11-13 19:09
dongyunfeng
Rank: 1
等 级:新手上路
帖 子:18
专家分:0
注 册:2006-11-14
收藏
得分:0 
我是菜鸟,求教大家,如何实现将VB中的数值赋值给MATLAB中的函数?我的联系方式
QQ:409104373 E-mail:dongyunfeng2004@126.com
2006-11-14 12:36
hitzhang
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:版主
威 望:21
帖 子:369
专家分:52
注 册:2006-9-24
收藏
得分:0 

不知道算得对不对


function dy=fy(x,y,z,C1,C2)
Gn=8.4e-13;n0=1.4e+24;tin=8e-12;k=0.015;a=6;
dy=(1/2)*a*Gn*(z-n0)-(k/tin)*(C1/x)*sin(y-C2);


function dx=fx(x,y,z,C1,C2)
Gn=8.4e-13;nth=2.0e+24;k=0.015;tin=8e-12;
dx=(1/2)*Gn*(z-nth)*x+(k/tin)*C1*cos(y-C2);



function dz=fz(x,z)
Gn=8.4e-13;e=1.6e-19;V=1.0e-16;ts=2e-9;Ith=54e-3;n0=1.4e+24;
dz=(1.3*Ith)/(e*V)-(1/ts)*z-Gn*(z-n0)*x^2;



function [T,X,Y,Z]=qwe()




C1=[8.2e-011, -1.4679e-007, 2.2247e-007, 2.3087e-007, -1.0047e-007, 3.153e-008, 3.2047e-007, 1.0044e-007, -1.7736e-007, 1.8828e-007, -1.9584e-007, -1.1957e-007, 2.6189e-007, 2.4234e-008, 1.6046e-007, -1.2331e-007, -1.9984e-007, 4.6585e-007, -2.9163e-007, 1.4481e-008, -1.2686e-007, 6.6135e-007, -1.1262e-007, 1.3977e-007, -6.5892e-008, -3.6625e-007, -2.467e-007, 7.4485e-008, 1.475e-007, 2.7209e-007, 2.4435e-007, 1.0191e-007, 8.4376e-008, 1.5446e-008, -1.4789e-007, -1.9714e-007, 3.0642e-007, -3.1694e-007, 1.3063e-007, 1.2397e-007, -5.0877e-008, 8.2291e-008, -4.9396e-008, 4.4198e-008, -6.6902e-010, -2.0217e-008, 1.4079e-008, 2.1972e-008, 7.8652e-009, 1.0089e-008, 2.0888e-008, -5.0502e-008, 1.75e-007, 8.8072e-007, 1.6731e-006, 2.1356e-006, 2.0616e-006, 1.4339e-006, 2.6416e-007, -1.1631e-006, -2.3767e-006, -2.8328e-006, -1.3912e-006, 3.3939e-006, 1.7953e-005, 0.00010396, 0.00069155, 0.0052414, 0.045308, 0.44673, 5.0174, 64.156, 933.44, 15451, 2.9058e+005, 6.2072e+006, 1.5046e+008, 4.1381e+009, 1.2671e+011, 4.4426e+011, 8.9674e+010, 1.8253e+010, 4.2933e+009, 1.1719e+009, 3.7112e+008, 1.3615e+008, 5.7867e+007, 2.8467e+007, 1.6195e+007, 1.0647e+007, 8.0882e+006, 7.0847e+006, 7.1595e+006, 8.3524e+006, 1.1181e+007, 1.7265e+007, 3.0646e+007, 6.2556e+007, 1.4634e+008, 3.9367e+008, 1.2152e+009, 4.2877e+009, 1.7323e+010, 7.8254e+010, 2.8099e+011, 2.251e+011, 8.7931e+010, 3.4687e+010, 1.5524e+010, 8.0087e+009, 4.7717e+009, 3.2836e+009, 2.609e+009, 2.391e+009, 2.5263e+009, 3.0734e+009, 4.3131e+009, 6.934e+009, 1.2876e+010, 2.7086e+010, 6.4164e+010,
];
C2=[0, -5.3281, -9.5561, -12.689, -14.733, -15.694, -15.576, -14.385, -12.126, -8.8054, -4.4275, 1.0022, 7.4783, 14.996, 23.549, 33.133, 43.744, 55.374, 68.021, 81.677, 96.34, 112, 128.66, 146.31, 164.94, 184.56, 205.15, 226.71, 249.24, 272.73, 297.17, 322.57, 348.91, 376.2, 404.42, 433.58, 463.66, 494.67, 526.59, 559.44, 593.19, 627.84, 663.4, 699.85, 737.2, 775.43, 814.54, 854.54, 895.4, 937.14, 979.74, 1023.2, 1067.5, 1112.7, 1158.7, 1205.6, 1253.3, 1301.8, 1351.2, 1401.4, 1452.4, 1504.3, 1556.9, 1610.4, 1664.7, 1719.7, 1775.6, 1832.3, 1889.8, 1948, 2007, 2066.8, 2127.4, 2188.8, 2250.9, 2313.8, 2377.5, 2441.9, 2506.9, 2558.9, 2593.8, 2628.8, 2664.6, 2701.4, 2739, 2777.5, 2816.9, 2857.1, 2898.3, 2940.3, 2983.2, 3026.9, 3071.5, 3116.9, 3163.2, 3210.3, 3258.3, 3307.1, 3356.7, 3407.1, 3458.4, 3510.5, 3563.4, 3617, 3669.2, 3712.4, 3751.3, 3790.2, 3829.9, 3870.5, 3911.9, 3954.2, 3997.3, 4041.3, 4086.2, 4131.9, 4178.4, 4225.8, 4274, 4323, 4372.7,
];

a=0.000000012;
b=0.000000024;
M=3000;
h=(b-a)/M;
t=linspace(a,b,3001);
x=zeros(1,M+1);
y=zeros(1,M+1);
z=zeros(1,M+1);
t=a:h:b;
x(1)=6.4164e+010;
y(1)=4372.7;
z(1)=2.2129e+024;
for j=1:M

kx1=h*feval('fx',x(j),y(j),z(j),C1(j),C2(j));
ky1=h*feval('fy',x(j),y(j),z(j),C1(j),C2(j));
kz1=h*feval('fz',x(j),y(j));


kx2=h*feval('fx',x(j)+kx1/2,y(j)+ky1/2,z(j)+kz1/2,C1(j),C2(j));
ky2=h*feval('fy',x(j)+kx1/2,y(j)+ky1/2,z(j)+kz1/2,C1(j),C2(j));
kz2=h*feval('fz',x(j)+kx1/2,z(j)+kz1/2);


kx3=h*feval('fx',x(j)+kx2/2,y(j)+ky2/2,z(j)+kz2/2,C1(j),C2(j));
ky3=h*feval('fy',x(j)+kx2/2,y(j)+ky2/2,z(j)+kz2/2,C1(j),C2(j));
kz3=h*feval('fz',x(j)+kx2/2,z(j)+kz2/2);


kx4=h*feval('fx',x(j)+kx3,y(j)+ky3,z(j)+kz3,C1(j),C2(j));
ky4=h*feval('fy',x(j)+kx3,y(j)+ky3,z(j)+kz3,C1(j),C2(j));
kz4=h*feval('fz',x(j)+kx3,z(j)+kz3);


x(j+1)=x(j)+(kx1+2*kx2+2*kx3+kx4)/6;
y(j+1)=y(j)+(ky1+2*ky2+2*ky3+ky4)/6;
z(j+1)=z(j)+(kz1+2*kz2+2*kz3+kz4)/6;
C1=[C1 x(j+1)];
C2=[C2 y(j+1)];

end
T=t;
X=x;
Y=y;
Z=z;

plot(t,x);

2006-11-15 15:29
快速回复:求助,matlab求解微分方程组
数据加载中...
 
   



关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.023507 second(s), 8 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved