求教：用matlab求解非线性微分方程组
本帖最后由 neverquit1988 于 2012-10-8 19:24 编辑大牛们：小弟想用matlab表示出下述关于蠕变非线性微分方程组的ε-t曲线，请各位大牛指教。如果有还需要提供的数据，请您提出。谢谢！ε, Η, ω, Φ,都是关于t的函数，初值是ε（0）=0， H(0)=0，Φ(0)=0，ω(0)=0.
把我命名的函数贴出来，避免重复劳动.function dy=creep(t,y)
dy=zeros(4,1)
a=%参数的值
s=110 %施加的应力
dy(1)=a(1)*sinh(a(2)*110*(1-y(2))/(1-y(3))/(1-y(4)))
dy(2)=a(4)*(a(1)*sinh(a(2)*110*(1-y(2))/(1-y(3))/(1-y(4))))/110/(1-(y(2)/a(5)))
dy(3)=a(6)/3/(1-dy(3))^4
dy(4)=a(3)*(a(1)*sinh(a(2)*110*(1-y(2))/(1-y(3))/(1-y(4))))
end
本帖最后由 neverquit1988 于 2012-9-24 08:50 编辑
请各位老师给予指导，谢谢! 这是个隐函数的微分方程组，我试着用matlab的ODE解，没得到理想的结果。 用1stOpt试了下，不知对否，参考下：
Constant a=;
Variable x=, y(4)=0;
Plot y1, y2, y3, y4;
ODEFunction y1'=a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4));
y2'=a4*(a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4)))/110/(1-(y2/a5));
y3'=a6/3/(1-y3')^4;
y4'=a3*(a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4)));
x y1(x) y2(x) y3(x) y4(x)
0 0 0 0 0
0.5 2.25094783175176E-5 0.00454556058462381 3.99655361132745E-6 0.000126512321456628
1 4.37986848632223E-5 0.00888905538320221 7.99310722265491E-6 0.000246166224762361
1.5 6.39821025524339E-5 0.0130479862547073 1.19896608339824E-5 0.000359605149946325
2 8.31592195833342E-5 0.0170376722967025 1.59862144453098E-5 0.000467388260696455
2.5 0.000101417009030962 0.0208715985338732 1.99827680566373E-5 0.000570004380675038
3 0.000118831919469203 0.0245616976334118 2.39793216679647E-5 0.00066788318161493
3.5 0.000135471487484506 0.0281185795153204 2.79758752792922E-5 0.00076140424629519
4 0.000151395654765219 0.0315517200117585 3.19724288906196E-5 0.000850904471112879
4.5 0.000166657852154546 0.0348696170394568 3.59689825019471E-5 0.000936684158896686
5 0.000181305898249052 0.0380799207759569 3.99655361132745E-5 0.00101901206939195
5.5 0.000195382749194802 0.0411895428660426 4.3962089724602E-5 0.00109812963341651
6 0.000208927128175809 0.0442047485861597 4.79586433359294E-5 0.001174254490839
6.5 0.000221974056938811 0.0471312350621782 5.19551969472569E-5 0.00124758347796182
7 0.000234555307016513 0.0499741979991265 5.59517505585844E-5 0.00131829516357728
7.5 0.000246699784715787 0.052738388890234 5.99483041699118E-5 0.00138655201275613
8 0.000258433861152835 0.0554281642904172 6.39448577812392E-5 0.00145250224177789
8.5 0.000269781656444062 0.0580475284397156 6.79414113925667E-5 0.00151628141539785
9 0.000280765285452746 0.0606001702855926 7.19379650038942E-5 0.00157801382804224
9.5 0.000291405071138886 0.0630894957649132 7.59345186152216E-5 0.00163781370292015
10 0.00030171973048161 0.0655186560559034 7.99310722265491E-5 0.00169578623698225
本帖最后由 neverquit1988 于 2012-10-16 09:06 编辑
shamohu 发表于 2012-9-28 16:28 static/image/common/back.gif
用1stOpt试了下，不知对否，参考下：x y1(x) y2(x) y3(x) y4(x)
0 0 0 0 0
0.5 2.25094783175176E-5 0.0045 ...
版主，您好！首先非常感谢您的帮助！您给的x的范围有点小，x相当于蠕变的时间在这里应该在（0,30000）之间，得出的x-y的关系应尽量满足这样一组数据：
x y1
3179.056255 0.22579144
10410.7112 0.58132454
14168.1403 0.839411788
16770.1533 1.000797756
18792.46387 1.16199758
20236.40159 1.290901597
21100.63687 1.419619471
22543.245 1.580633152
23981.86435 1.837975828
25133.29167 2.031006103
25994.86777 2.223943306
26853.78469 2.481099839
27426.83916 2.641834306
28285.75607 2.898990839
28563.64096 3.188070894
28844.18503 3.412931619
29124.7291 3.637792344
29402.61398 3.926872399
29683.15805 4.151733124
29671.19171 4.440720107
29661.88456 4.665487761
29651.24782 4.92236508
29930.4623 5.17933547
29627.31515 5.500339046
29610.03045 5.917764689
版主，您好，得到的x-y曲线应该与这个曲线很接近的。谢谢您！ 版主很热心，很感谢。 dbx12358 发表于 2012-9-28 19:30 static/image/common/back.gif
版主很热心，很感谢。
希望版主能够继续提供帮助，1stopt我只有网上1.5版本的，被阉割了，不能解微分方程组。 本帖最后由 feiyuzhen 于 2012-10-5 09:51 编辑
用maple算的，代码如下 restart;
sys := diff(epsilon(t), t) = A*sinh(B*sigma*(1-H(t))/((1-Phi(t))*(1-omega(t)))), diff(H(t), t) = h*(diff(epsilon(t), t))*(1-H(t)/H1)/sigma, diff(Phi(t), t) = (1/3)*Kc*(1-Phi(t))^4, diff(omega(t), t) = C*(diff(epsilon(t), t));
A := 0.22959353e-9; B := .11734557; C := 5.6204022; h := 22097.986; H1 := .43749936; Kc := 0.23978555e-4; sigma := 110;
IC := epsilon(0) = 0, H(0) = 0, Phi(0) = 0, omega(0) = 0;
dsolve({IC, sys}, numeric, output = listprocedure);
plot(`~`([%, %, %, %]), 0 .. 30000, thickness = 2, gridlines = true, legend = );
feiyuzhen 发表于 2012-10-5 09:47 static/image/common/back.gif
用maple算的，代码如下
谢谢版主的热心帮助，这两天放假回老家了，不好意思。我这两天研究一下，呵呵! 本帖最后由 neverquit1988 于 2012-10-8 15:48 编辑
feiyuzhen 发表于 2012-10-5 09:47 static/image/common/back.gif
用maple算的，代码如下
版主是否能提供用MATLAB的解决办法，我电脑上没有maple。我用matlab的ODE15i计算老是出错，代码如下：a=%参数的值
s=110 %施加的应力
odefun=@(t,y,dy)[dy(1)-a(1)*sinh(a(2)*110*(1-y(2))/(1-y(3))/(1-y(4)))
dy(2)-a(4)*(a(1)*sinh(a(2)*110*(1-y(2))/(1-y(3))/(1-y(4))))/110/(1-(y(2)/a(5)))
dy(3)-a(6)/3/(1-dy(3))^4
dy(4)-a(3)*(a(1)*sinh(a(2)*110*(1-y(2))/(1-y(3))/(1-y(4))))]
t0=0;%自变量初值
y0=';
fix_y0=ones(4,1);
dy0=';
fix_dy0=zeros(4,1);
=decic(odefun,t0,y0,fix_y0,dy0,fix_dy0)
=ode15i(odefun,,y02,dy02)
网上maple有很多，你下一个就行，matlab有段时间不用了，你看看薛定宇的那本书应该可以解决 feiyuzhen 发表于 2012-10-9 09:58 static/image/common/back.gif
网上maple有很多，你下一个就行，matlab有段时间不用了，你看看薛定宇的那本书应该可以解决 ...
好的，谢谢了！有些结果是能算出来，但结果不合适，还是谢谢了！我怎么给你们金币呢？找了好久也没找到如何支付。 本帖最后由 zzz700 于 2012-10-10 11:12 编辑
neverquit1988 发表于 2012-9-28 19:17 static/image/common/back.gif
版主，您好！首先非常感谢您的帮助！您给的x的范围有点小，x相当于蠕变的时间在这里应该在（0,30000）之间 ...
当微分方程组的各项参数确定之后，微分方程组的解也就确定了。
如：
shamohu 发表于 2012-9-28 16:28 static/image/common/back.gif
用1stOpt试了下，不知对否，参考下：x y1(x) y2(x) y3(x) y4(x)
0 0 0 0 0
0.5 2.25094783175176E-5 0.0045 ...
而你称“得出的x-y的关系应尽量满足这样一组数据”
是不是指调整a、s的值，
（a=%参数的值
s=110 %施加的应力），
使微分方程组得出的x-y的关系应尽量满足这样一组数据？
本帖最后由 neverquit1988 于 2012-10-12 17:11 编辑
zzz700 发表于 2012-10-10 11:10 static/image/common/back.gif
当微分方程组的各项参数确定之后，微分方程组的解也就确定了。
如：
不是，是（x，y1）的关系应该尽量趋近与我给的数据。a和s的值不需要调整都是给定的，就是解一下这个微分方程组得出x，y的曲线。根据4楼版主的回复，我把方程改了一下。您如果有1stopt，麻烦您帮我运行一下这个代码，求得x，y1的对应数值。顺带问您一下，你的1stopt是正版的吗？Constant a=;
Variable x=, y(4)=0;
Plot y1, y2, y3, y4;
ODEFunction y1'=a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4));
y2'=a4*(a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4)))/110/(1-(y2/a5));
y3'=a6/3/(1-y3)^4;
y4'=a3*(a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4)));
本帖最后由 zzz700 于 2012-10-12 19:49 编辑
neverquit1988 发表于 2012-10-11 07:47 static/image/common/back.gif
不是，是（x，y1）的关系应该尽量趋近与我给的数据。a和s的值不需要调整都是给定的，就是解一下这个微分方 ...
ODEFunction y1'=a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4));
y2'=a4*(a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4)))/110/(1-(y2/a5));
y3'=a6/3/(1-y3)^4;
y4'=a3*(a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4)));
matlab的ode45
ODEFunction y1'=a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4));
y2'=a4*(a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4)))/110/(1-(y2/a5));
y3'=a6/3/(1-y3')^4;
y4'=a3*(a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4)));
zzz700 发表于 2012-10-12 19:53 static/image/common/back.gif
ODEFunction y1'=a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4));
y2'=a4*(a1*sinh(a2*110*(1-y2)/(1- ...
这个是用1stopt的结果？ zzz700 发表于 2012-10-12 19:47 static/image/common/back.gif
ODEFunction y1'=a1*sinh(a2*110*(1-y2)/(1-y3)/(1-y4));
y2'=a4*(a1*sinh(a2*110*(1-y2)/(1- ...
麻烦您把matlab运行的代码贴出来看看，我也用ode45解，为什么没解出你的结果？
页:
[1]