一个求解常微分方程组的程序老是出错，求各位大神指导下，谢谢了

chourhy

方程是这样的：dx/dt=cosa   dy/dt=sina  da/dt=m  dm/dt=-10sina
边界条件是：   t=0时，x=0,y=0,m=0
                        t=1时，y=0,m=0
本人运行的程序如下：
Clear[pi,gm,m,c]
pi = N[\[Pi]];
pe=10;
gm = 0.53 ;
de[y3_,y4_]:={y1'[t]\[Equal]Cos[y3[t]],y2'[t]\[Equal]Sin[y3[t]],
   y3'[t]\[Equal]y4[t],y4'[t]\[Equal]-pe*Sin[y3[t]]}
leftBC[m_]:={y1[0]\[Equal]0,y2[0]\[Equal]0,y3[0]\[Equal]m,y4[0]\[Equal]0}
soln:=NDSolve[Flatten[Append[de[y3,y4],leftBC[m]]],{y1,y2,y3,y4},{t,0,1},
   MaxSteps\[Rule]2000]
endpt[m_]:={y1[t],y2[t],y3[t],y4[t]}/.First[
       NDSolve[Flatten[Append[de[y3,y4],leftBC[m]]],{y1[t],y2[t],y3[t],
         y4[t]},{t,0,1},MaxSteps\[Rule]2000]]/.t\[Rule]1;
Clear[m]
rts:=FindRoot[{endpt[m][[2]]\[Equal]0},{m,{gm,0.99gm}},AccuracyGoal\[Rule]6,
   MaxIterations\[Rule]1000]
endpt[m/.rts]
m=m/.rts
{yy1[t_],yy2[t_],yy3[t_],yy4[t_]}={y1[t],y2[t],y3[t],y4[t]}/.First[soln];
ParametricPlot[Evaluate[{yy1[t],yy2[t]}/.soln/.rts],{t,0,1},
PlotRange\[Rule]All,AspectRatio\[Rule]Automatic,PlotPoints\[Rule]100]
但是老是弹出错误提示：NDSolve::ndinnt: Initial condition m is not a number or a rectangular array of numbers
实在不清楚哪里有问题，所以来求助各位高手，希望不吝赐教  
谢谢了

xhn

endpt函数只能接受数值变量，Findroot会对方程进行符号代入化简，导致NDSolve以符号m做为初值，计算出错。
Findroot好像不能以这种方式给定初值。
程序按如下方式修改，不报错，能得到结果。
对不对我就不确定了，说实话没看懂你的程序……

endpt[m_?NumericQ] := ({y1[t], y2[t], y3[t], y4[t]} /.
      First[NDSolve[
        Flatten[Append[de[y3, y4], leftBC[m]]], {y1[t], y2[t], y3[t],
         y4[t]}, {t, 0, 1}, MaxSteps -> 2000]] /. t -> 1)[[2]];
Clear[m]
rts := FindRoot[{endpt[m] == 0}, {m, gm, 0.99 gm}, AccuracyGoal -> 6,
  MaxIterations -> 1000]

chourhy

回复 2 # xhn 的帖子

非常感谢你的解答
但是我按照你的意思改了过来为什么还是出错啊。。。
错误提示是：FindRoot::nlnum: "The function value {endpt[{0.53,0.5247}]} is not a list of numbers with dimensions {1} at {m} = {{0.53,0.5247}}"
不知道是哪里出了问题  
能解答下吗  谢谢了

chourhy

回复 2 # xhn 的帖子

早上起来又试了一下 可以算了  谢谢你了