求助，大难题，微分方程组振荡

debuger

这是啥原因？刚性问题？还是稳定性问题，谢谢指教。其中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);

用定步长：[t,x]=ode45(@fun2,tspan,[30,250000,2.55,343,30/sqrt(1.4*287*343)])



用变步长：[t,x]=ode45(@fun2,[0 2.41],[30,250000,2.55,343,30/sqrt(1.4*287*343)])

donkeyxu

??? function xdot = fun2(t,x)
|
Error: Function definitions are not permitted at the prompt or in scripts

debuger

<P>我是MATLAB7.0,是是可以通过的啊</P>

happy

不知道你想问什么？我这里计算结果不管是变步长还是定步长都一样的
都是你给出来的第一个图

debuger

对啊，后部强烈振荡是为什么？谢谢这是一个异形喷管的的流道，后部强烈扩展，是应该有影响的

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(:);

happy

呵呵，这个就不是matlab的事情了，具体问题就要你分析了，matlab不是万能

1. 分析一下无理模型是否有可能出现这种情况，如果可能那就没问题了
2. 如果不可能检查你的数学模型是否有正确，正确的话他是否能放映物理模型
3. 数学模型没有问题的话，那就检查你的function吧，是否什么地方弄错了


流场之类的我可看不懂

happy

另外你后面又加了个程序是什么意思？

debuger

global da;
global a;我说这个程序是可以有我开始说的那个问题的，谢谢

debuger

就是，定步长的结果最后扩展部分混乱，<BR>变步长的，整个解都混乱。应该是数值方法的问题啊，我认为。。。

happy

用 test运行结果还是一样

siyanger

我也算了一下，也都是第一个图啊。
变步长的，整个解都混乱，怎么混乱呢？第二个图，以振动方程来看，结果还是很有规律的呢。我求解的振动微分方程就希望出现这样的结果呢！

debuger

自己UP

debuger

<P>自己顶上去</P>