Клен. Dsolve и функции

Question

У меня есть дифференциальное уравнение, которое решает моя программа:

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

После решения мы имеем следующее:

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

мне нужно использовать этот альфа (Т) и бета (т) в следующим образом:

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. Те более эффективны для вызова много раз, в отличие от вашего оригинала, который формировал последовательность значений для каждого числового значения t, а затем приходилось выбирать определенный операнд. Он также более устойчив, поскольку он не чувствителен к позициям (который может измениться из-за нового лексикографического упорядочения, если вы изменили имена ваших переменных).

Вышеупомянутое заняло около 16 секунд на Intel i7. Это не быстро. Одна из причин заключается в том, что числовые интегралы вычисляются с большей точностью, чем это необходимо для построения графика. Итак, давайте перезагрузим (чтобы обеспечить справедливое время) и перекомпонуем с ослабленными допусками на числовые интегрирования, выполненные для 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 не являются функциями, для которых у нас есть явные функции - - их точность тесно связана с численным решением. Это не совсем то же самое, что упрощенно использовать числовую процедуру 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;