0

我的程序解决了一个微分方程：

p := dsolve({ic, sys}, numeric, method = rosenbrock);

解决后我们有以下内容：

print(p(10));
      [t = 10., alpha(t) = HFloat(0.031724302221312055), beta(t) = HFloat(0.00223975915581258)]

我需要使用这个 alpha(t) 和 beta(t)，如下所示：

a := t->exp( int(alpha(t)),x=0..t) );
b := t->exp( int(beta(t)),x=0..t) )

并绘制一个情节：

odeplot(p, [[t, a(t)], [t, b(t)]], 0 .. 20, thickness = 2, numpoints = 500, color = [red, blue])

这样做的第一件事是：

p := dsolve({sys, ic},  numeric, method=rosenbrock); 
alpha := t->rhs(p(t)[2] );
beta := t->rhs(p(t)[3;
a := t->exp( int(alphat)),x=0..t) );
b := t->exp( int(betat)),x=0..t) );
odeplot(p, [[t, a(t)], [t, b(t)]], 0 .. 20, thickness = 2, numpoints = 500, color = [red, blue])

但是代码不起作用，是的，显然，应该采取不同的行动。

4

1 回答 1

1

首先，让我们尝试通过一些语法和用法更改来让您的方法发挥作用。

您没有提供示例的详细信息，因此我构建了一个微分方程系统sys和初始条件ic，以便可以执行您的计算命令。

restart:

sys := diff(alpha(t),t) = -1/200*beta(t),
       diff(beta(t),t) = -1/200*alpha(t) - 1/Pi^2:
ic := alpha(0)=0, beta(0)=-1:

p := dsolve({ic, sys}, numeric, method = rosenbrock, output=listprocedure):

alphat := eval(alpha(t),p):
betat := eval(beta(t),p):

a := unapply( exp( Int(alphat , 0..t) ), t, numeric):
b := unapply( exp( Int(betat , 0..t) ), t, numeric):

evalf(a(20.0)), evalf(b(20.0));

                                              -18
                   5.347592595, 3.102016550 10   

st := time():
P := plots:-odeplot(p, [[t, a(t)], [t, b(t)]], 0 .. 20,
               thickness = 2, numpoints = 50, color = [red, blue]):
( time() - st )*`seconds`;

                           16.770 seconds

P;

在此处输入图像描述

我使用output=listprocedure这样我可以将正确的程序从解决方案分配p到alphat和betat。与您的原始数据相比，这些数据多次调用效率更高，原始数据为每个数值形成一系列t值，然后必须选择某个操作数。它也更加健壮，因为它对位置不敏感（如果您更改了变量的名称，这可能会由于新的字典顺序而改变）。

以上在 Intel i7 上花费了大约 16 秒。这并不快。原因之一是数值积分的计算精度高于绘图所需的精度。a因此，让我们重新开始（以确保公平的时间）并重新计算和的数值积分的宽松容差b。

restart:
sys := diff(alpha(t),t) = -1/200*beta(t),
       diff(beta(t),t) = -1/200*alpha(t) - 1/Pi^2:
ic := alpha(0)=0, beta(0)=-1:
p := dsolve({ic, sys}, numeric, method = rosenbrock, output=listprocedure):
alphat := eval(alpha(t),p):
betat := eval(beta(t),p):

a := unapply( exp( Int(alphat , 0..t, epsilon=1e-5) ), t, numeric):
b := unapply( exp( Int(betat , 0..t, epsilon=1e-5) ), t, numeric):

evalf(a(20.0)), evalf(b(20.0));

                                              -18
                   5.347592681, 3.102018090 10   

st := time():
P := plots:-odeplot(p, [[t, a(t)], [t, b(t)]], 0 .. 20,
               thickness = 2, numpoints = 50, color = [red, blue]):
( time() - st )*`seconds`;

                            0.921 seconds

您可以检查这是否给出了相同的图。

现在让我们扩充这个例子，这样数字积分就可以dsolve,numeric自己计算了。我们可以通过使用微积分来做到这一点。通过这种方式，我们将其留给数字 ode 求解器进行自己的误差估计、步长控制、刚度或奇异性检测等。

请注意，被积函数alphat和betat不是我们有明确函数的函数——它们的准确性与数值 ode 求解密切相关。这与简单地使用数值 ode 例程替换数值求积例程来解决我们期望直接计算到任何所需精度（包括任何奇点的任一侧）的被积函数的问题并不完全相同。

restart:

sys := diff(alpha(t),t) = -1/200*beta(t),
       diff(beta(t),t) = -1/200*alpha(t) - 1/Pi^2,
       diff(g(t),t) = alpha(t), diff(h(t),t) = beta(t):
ic := alpha(0)=0, beta(0)=-1,
      g(0)=0, h(0)=0:

p := dsolve({ic, sys}, numeric, method = rosenbrock, output=listprocedure):

alphat := eval(alpha(t),p):
betat := eval(beta(t),p):
gt := eval(g(t),p):
ht := eval(h(t),p):

exp(gt(20.0)), exp(ht(20.0));

                                                   -18
              5.34759070530497, 3.10201330730572 10   

st := time():
P := plots:-odeplot(p, [[t, exp(g(t))], [t, exp(h(t))]], 0 .. 20,
               thickness = 2, numpoints = 50, color = [red, blue]):
( time() - st )*`seconds`;

                            0.031 seconds
P;

在此处输入图像描述

于 2013-07-21T04:33:05.687 回答