楓葉。 Dsolve和功能

Question

我有我的程序解微分方程：

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

解決，我們有以下後：

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

我需要用這個α（t）和β（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]) 

但代碼不能正常工作，並且是的，很明顯，應該採取不同的。

Answer 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。那些調用多次的效率更高，而不是你的原始值，它爲每個數值形成了一系列值，然後必須選取某個操作數。它也更強大，因爲它對位置不敏感（如果您更改了變量的名稱，可能會由於新的字典順序而改變）。

以上在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本身計算。我們可以通過使用微積分來做到這一點。通過這種方式，我們將它留給數字頌求解器來做自己的誤差估計，步長控制，剛度或奇點檢測等。

請注意，被積函件alphat和betat不是我們有明確函數的函數 - - 它們的準確性與數字解析密切相關。這與簡單地使用數字代碼例程替換數字正交例程的結果並不完全一樣，因爲我們希望能夠直接計算任意期望的精度（包括任何奇點的任一側）的被積函數問題。

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;