Не удалось решить простой ODE с Octave

Я новичок в Octave, поэтому я пытаюсь заставить работать несколько простых примеров, прежде чем переходить к более сложным проектам.

Я пытаюсь решить ODE dy/dx = a*x+b, но безуспешно. Вот код:

%Funzione retta y = a*x + b. Ingressi: vettore valori t; coefficienti a,b
clear all;
%Inizializza argomenti
b = 1;
a = 1;
x = ones(1,20);
function y = retta(a, x, b) %Definisce funzione
y = ones(1,20);
y = a .* x .+ b;
endfunction
%Calcola retta
x = [-10:10];
a = 2;
b = 2;
r = retta(a, x, b)
c = b;
p1 = (a/2)*x.^2+b.*x+c  %Sol. analitica di dy/dx = retta %
plot(x, r, x, p1);
% Risolve eq. differenziale dy/dx = retta %
y0 = b; x0 = 0;
p2 = lsode(@retta, y0, x)

И вывод:

retta3code
r =

 -18  -16  -14  -12  -10   -8   -6   -4   -2    0    2    4    6    8   10 12   14   16   18   20   22

p1 =

Columns 1 through 18:

82    65    50    37    26    17    10     5     2     1     2     5    10    17    26    37    50    65

Columns 19 through 21:

82   101   122

error: 'b' undefined near line 9 column 16
error: called from:
error:   retta at line 9, column 4
error: lsode: evaluation of user-supplied function failed
error: lsode: inconsistent sizes for state and derivative vectors
error:   /home/fabio/octave_file/retta3code.m at line 21, column 4

Итак, функция retta работает правильно в первый раз, но не работает при использовании в lsode, Почему это происходит? Что нужно изменить, чтобы код работал?

Источник

2 ответа

Решение

Каким-то образом вы все еще пропускаете некоторые важные части истории. Чтобы решить ODE y'=f(y,x) вам нужно определить функцию

function ydot = f(y,x)

где ydot имеет те же размеры, что и y, оба должны быть векторами, даже если они имеют размерность 1. x это скаляр По какой-то традиционной причине, lsode (код FORTRAN, используемый в нескольких пакетах решателей), предпочитает менее используемый порядок (y,x), в большинстве учебников и других решателей вы найдете порядок (x,y),

Затем получить образцы раствора ylist по точкам выборки xlist ты звонишь

ylist = lsode("f", y0, xlist)

где xlist(1) это начальное время.

Внутренние органы f не зависят от списка образцов list и какой у него размер. Это отдельная проблема, которую вы можете использовать мульти-оценки, чтобы вычислить точное решение с чем-то вроде

yexact = solexact(xlist)

Для передачи параметров используйте анонимные функции, как в

function ydot = f(y,x,a,b)
    ydot = [ a*x+b ]
end

a_val = ...
b_val = ...
lsode(@(y,x) f(y,x,a_val, b_val), y0, xlist)
Источник

Код, измененный ниже, работает, но я бы предпочел иметь возможность определять параметры a а также b из функции, а затем передать их rdot в качестве аргументов.

x = [-10,10];
a = 1;
b = 0;
c = b;
p1 = (a/2).*(x.^2)+b.*x+c  %Sol. analitica di dy/dx = retta %
function ydot = rdot(ydot, x)
a = 1;
b = 0;
ydot = ones(1,21);
ydot = a.*x .+ b;
endfunction
y0 = p1(1); x0 = 0;
p2 = lsode("rdot", y0, x, x0)'
plot(x, p1, "-k", x, p2, ".r");
Источник
Другие вопросы по тегам