feddy ha scritto:Visto che lo shooting ist klar...
Prima cosa: ridurre la ODE ad un sistema di ODEs del prim'ordine.
Poi, $u_n \approx u(t_n)$. E ovviamente con $t_{n+1}$ si intende $t_n+h$, ossia l'istante temporale successivo al tempo $n$.
La "mesh" spaziale qui sarà un semplice "linspace" (visto che il metodo è supposto essere a passo costante) tra $0$ e $L$.
Ora, i metodi Runge-Kutta sono facili da implementare (è un ciclo for), l'unica cosa che **potrebbe** complicare la situazioni è se il metodo è implicito (o semiimplicito): per capirlo, guarda il tableau che hai davanti (e magari apri un testo introduttivo). Chiaramente non è implicito, sarà semiimplicito oppure esplicito. La differenza in un caso oppure nell'altro potrebbe essere che ad ogni passo è richiesta la risoluzione di un sistema non lineare per trovare i nuovi $K_i$, **ma**, visto che la ODE è lineare, anche nel caso di un metodo semiimplicito avrai a che fare con sistemi lineari. In particolare, puoi fattorizzare la matrice di iterazione del metodo una volta per tutte in LU (ad esempio) e risolvere tutti i sistemi con la matrice fattorizzata. Questa cosa può sembrarti in un certo senso "frivola", ma non lo è a mio avviso, visto che con un metodo di shooting ti è richiesto di risolvere più volte lo stesso problema.
Qual è l'ordine di tale metodo? E'
A-Stabile?
Per verificare l'implementazione puoi risolverla anche con differenze finite centrate, visto che il problema è molto semplice e ti basta solamente discretizzare la derivata seconda come al solito e il problema risulta (chiamo l'incognita $\varphi =u$ per comodità)
\[ -k \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}+b u_i = 1\]
con $u_1=0$ e la condizione di Neumann al bordo da imporre ad esempio con un nodo fantasma (o anche differenze finite decentrate).
Ok, sono riuscito ad implementare il metodo di Runge-Kutta utilizzando l'array dato dalla traccia.
Il problema adesso è lo shooting; o meglio, pensavo non fosse un problema ma lo è diventato adesso perché prima ero solito risolvere i problemi con lo shooting attraverso il tool ode45 di MatLab... E chiaramente suppongo che adesso ci voglia un altro approccio.
Cioè mi spiego. Questa è la risoluzione che avrei adottato col tool ode45.
Anzitutto creo una function che generi il mio sistema
Testo nascosto, fai click qui per vederlo
- Codice:
function dy = sistema(t,y)
global b k
dy = y(:);
dy(1) = y(2);
dy(2) = (b/k)*y(1)-1/k;
end
Poi creo una function che risolve questo sistema di ODE con il parametro di shooting "s" che viene fuori dal metodo di bisezione.
Testo nascosto, fai click qui per vederlo
- Codice:
function f = FunConv (s)
L = 1;
N = 100;
xspan = linspace(0,L,N);
options = odeset ('RelTol', 1e-10, 'AbsTol', 1e-10);
[t,y] = ode45(@sistema, xspan, [0 s], options);
f = y(end,2) + 1;
end
con questo che è, diciamo, il Main
Testo nascosto, fai click qui per vederlo
- Codice:
clc; clear;
global b k
b = 1;
k = 0.5;
N = 100;
L = 1;
xspan = linspace(0,L,N);
s_vect = linspace ( -5, 5, 100);
for is=1:length(s_vect)
s = s_vect(is);
f(is) = FunConv(s);
end
figure(1);
plot (s_vect,f, 'r.-', 'markersize', 8); grid on;
title ('Gli zeri della funzione y(end,2)+1');
xlabel('Parametro s');
ylabel('f = y(end,2)+1');
s1 = 0;
s2 = 2;
ys1 = FunConv(s1);
ys2 = FunConv(s2);
if ys1*ys2 > 0
error('Ridefinisci i parametri!');
end
ym = 1;
Tol = 1e-12;
while abs(ym) > Tol
xm = (s1+s2)/2;
ym = FunConv(xm);
if ym*ys1 < 0
s2 = xm;
else
s1 = xm;
end
end
disp(['Il valore di s cercato è ', num2str(xm), '.'])
options = odeset('RelTol', 1e-10, 'AbsTol', 1e-10);
[t,y] = ode45(@sistema, xspan, [0 xm], options);
figure(2);
plot (xspan, y(:,1), 'b'); grid on;
title('Soluzione');
xlabel('x');
ylabel('\phi(x)');
Cioè cosa faccio; prima creo un array di parametri "s" in modo tale da passarli alla function FunConv e plottare una funzione che è la funzione differenza tra la derivata prima nota (data dal problema) e quella calcolata col parametro di shooting "s".
Faccio tutto questo per inizializzare il ciclo di shooting, cioè faccio tutto questo per assegnare i due valori iniziali, per poter entrare nel metodo di bisezione in maniera corretta.
Dopodiché eseguo lo shooting, chiamando sempre la funzione FunConv e passandole di volta in volta il parametro s = xm finché non ho che la differenza tra derivata in x = L calcolata con shooting e derivata in x = L nota sia circa pari a zero.
Quindi plotto tutto.
Adesso quello che mi chiedo è posso utilizzare questo ragionamento di base anche per il metodo RK?
