Fehler beim Lösen einer einfachen ODE mit Octave

Question

Ich bin neu in Octave, also versuche ich ein paar einfache Beispiele zu erstellen, bevor ich zu komplexeren Projekten übergehe.

Ich versuche, die ODE dy/dx = a*x+b zu lösen, aber ohne Erfolg. Hier ist der Code:

%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) 

Und der Ausgang ist:

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 

So wird die Funktion retta richtig beim ersten Mal funktioniert, aber es funktioniert nicht, wenn sie in lsode verwendet. Warum passiert das? Was muss geändert werden, damit der Code funktioniert?

Answer 1

Irgendwie vermissen Sie noch einige wichtige Teile der Geschichte. Um eine ODE zu lösen y'=f(y,x) Sie eine Funktion

function ydot = f(y,x) 

wo ydot die gleichen Abmessungen wie y hat definieren müssen, haben beide Vektoren sein, auch f sie von Dimension 1. x ist ein Skalar. Aus irgendeinem traditionellen Grund bevorzugt lsode (ein FORTRAN-Code, der in mehreren Solver-Paketen verwendet wird) die weniger verwendete Reihenfolge (y,x), in den meisten Lehrbüchern und anderen Solvern finden Sie die Reihenfolge (x,y).

Dann Lösungsproben ylist über Abtastpunkte xlist Sie

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

erhalten nennen, wo xlist(1) die Anfangszeit ist.

Die Interna von f sind unabhängig von der Beispielliste list und welche Größe es hat. Es ist eine andere Frage, die Sie Multi-Auswertung können die exakte Lösung mit so etwas wie

yexact = solexact(xlist) 

---

Um pass parameters verwenden anonymous functions, wie in

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) 

Answer 2

Der unten geänderte Code funktioniert, aber ich möchte die Parameter a und b außerhalb der Funktion definieren und sie dann als Argumente an rdot übergeben.

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");