求助，大难题，微分方程组振荡
这是啥原因？刚性问题？还是稳定性问题，谢谢指教。其中da/a是一个 微分计算，不过a的函数a(z)已经知道，以上可变系数方法对不对？函数定义：
function xdot = fun2(t,x)
global da;
global a;
if (t>=0 )&(t<0.43379)
a=(-0.57735 * t + 0.810) ^ 2 * 3.14159265235 - (0.267949 * t + 0.240) ^ 2 * 3.14159265235 ;
da=2*-0.57735*(-0.57735 * t + 0.810) * 3.14159265235 - 2*(0.267949 * t + 0.240)*0.267949 * 3.14159265235 ;
end
if (t>=0.43379 )&(t<0.73826)
a=((-sqrt(0.420 ^ 2 - (t - 0.64379) ^ 2) + 0.92328) ^ 2 * 3.14159265235) - (0.267949 * t + 0.240) ^ 2 * 3.14159265235;
da=2*(-sqrt(0.420^2-(t-0.64379)^2)+0.92328)*(2*t-1.28758)/(2*sqrt(-t^2+1.287582*t-0.2380655841))- (0.267949 * t + 0.240)*0.267949 * 2 * 3.14159265235;
end
if (t>=0.73826 )&(t<2.090)
a=(0.230868 * (t - 2.410) + 0.900) ^ 2 * 3.14159265235 - (0.267949 * t + 0.240) ^ 2 * 3.14159265235;
da= 2*(0.230868 * (t - 2.410) + 0.900)*0.230868 * 3.14159265235 - 2*(0.267949 * t + 0.240)*0.267949 * 3.14159265235;
end
if (t>=2.090 )&(t<2.410)
a = (0.230868 * (t - 2.410) + 0.900) ^ 2 * 3.14159265235 - (-(t - 2.090) + 0.800) ^ 2 * 3.14159265235;
da=2*(0.230868 * (t - 2.410) + 0.900) *0.230868*3.14159265235 - 2*(-(t - 2.090) + 0.800) * 3.14159265235;
end
k=1.4;
m=x(1)/sqrt(k*287*x(4));
xdot(1) = x(1)*(1/(m^2-1))*(da/a);
xdot(2) = x(2)*(k*m^2/(1-m^2))*(da/a);
xdot(3) = x(3)*(k*m^2/(1-m^2))*(da/a);
xdot(4) = x(4)*((k-1)*m^2/(1-m^2))*(da/a);
xdot(5)=x(5)*(-2-(k-1)*m^2)/(2*(1-m^2))*(da/a);
xdot=xdot(:); % To make xdot a column
% End of FUN2.M
调用：
t0=0;tf=2.41;nsteps=1000;
tspan = linspace(t0, tf, nsteps);
用定步长：=ode45(@fun2,tspan,)
用变步长：=ode45(@fun2,,) ??? function xdot = fun2(t,x)
|
Error: Function definitions are not permitted at the prompt or in scripts <P>我是MATLAB7.0,是是可以通过的啊</P>
回复：（debuger）求助，大难题，微分方程组振荡
不知道你想问什么？我这里计算结果不管是变步长还是定步长都一样的都是你给出来的第一个图 对啊，后部强烈振荡是为什么？谢谢这是一个异形喷管的的流道，后部强烈扩展，是应该有影响的 以下这个是能跑的，有上述问题的
function xdot = test(t,x)
global da;
global a;
if (t>=0 )&(t<0.43379)
a=(-0.57735 * t + 0.810) ^ 2 * 3.14159265235 - (0.267949 * t + 0.240) ^ 2 * 3.14159265235 ;
da=2*-0.57735*(-0.57735 * t + 0.810) * 3.14159265235 - 2*(0.267949 * t + 0.240)*0.267949 * 3.14159265235 ;
end
if (t>=0.43379 )&(t<0.73826)
a=((-sqrt(0.420 ^ 2 - (t - 0.64379) ^ 2) + 0.92328) ^ 2 * 3.14159265235) - (0.267949 * t + 0.240) ^ 2 * 3.14159265235;
da=2*(-sqrt(0.420^2-(t-0.64379)^2)+0.92328)*(2*t-1.28758)/(2*sqrt(-t^2+1.287582*t-0.2380655841))- (0.267949 * t + 0.240)*0.267949 * 2 * 3.14159265235;
end
if (t>=0.73826 )&(t<2.090)
a=(0.230868 * (t - 2.410) + 0.900) ^ 2 * 3.14159265235 - (0.267949 * t + 0.240) ^ 2 * 3.14159265235;
da= 2*(0.230868 * (t - 2.410) + 0.900)*0.230868 * 3.14159265235 - 2*(0.267949 * t + 0.240)*0.267949 * 3.14159265235;
end
if (t>=2.090 )&(t<2.410)
a = (0.230868 * (t - 2.410) + 0.900) ^ 2 * 3.14159265235 - (-(t - 2.090) + 0.800) ^ 2 * 3.14159265235;
da=2*(0.230868 * (t - 2.410) + 0.900) *0.230868*3.14159265235 - 2*(-(t - 2.090) + 0.800) * 3.14159265235;
end
k=1.4;
m=x(1)/sqrt(k*287*x(4));
xdot(1) = x(1)*(1/(m^2-1))*(da/a);
xdot(2) = x(2)*(k*m^2/(1-m^2))*(da/a);
xdot(3) = x(3)*(k*m^2/(1-m^2))*(da/a);
xdot(4) = x(4)*((k-1)*m^2/(1-m^2))*(da/a);
xdot(5)=x(5)*(-2-(k-1)*m^2)/(2*(1-m^2))*(da/a);
xdot=xdot(:);
回复：（debuger）对啊，后部强烈振荡是为什么？谢谢这...
呵呵，这个就不是matlab的事情了，具体问题就要你分析了，matlab不是万能1. 分析一下无理模型是否有可能出现这种情况，如果可能那就没问题了
2. 如果不可能检查你的数学模型是否有正确，正确的话他是否能放映物理模型
3. 数学模型没有问题的话，那就检查你的function吧，是否什么地方弄错了
流场之类的我可看不懂
回复：（debuger）求助，大难题，微分方程组振荡
另外你后面又加了个程序是什么意思？ global da;global a;我说这个程序是可以有我开始说的那个问题的，谢谢 就是，定步长的结果最后扩展部分混乱，<BR>变步长的，整个解都混乱。应该是数值方法的问题啊，我认为。。。
回复：（debuger）global da;global a;我说这个程序是...
用 test运行结果还是一样回复：（happy）回复：（debuger）求助，大难题，微分...
我也算了一下，也都是第一个图啊。变步长的，整个解都混乱，怎么混乱呢？第二个图，以振动方程来看，结果还是很有规律的呢。我求解的振动微分方程就希望出现这样的结果呢！ 自己UP <P>自己顶上去</P>
页:
[1]