Sto risolvendo questa equazione di convezione lineare non stazionaria
$ (d\varphi)/dt = a(d\varphi)/dx $
nell'intervallo $ [0,L] $ con condizioni al contorno periodiche e condizione iniziale
$ \varphi(x,0) = (x-L)^2cos((2\pix)/L) $.
La traccia richiede di discretizzare la derivata spaziale in questo modo
$ 1/6\varphi_(i-1)^' + 2/3\varphi_i^' + 1/6\varphi_(i+1)^' = (\varphi_(i+1)-\varphi_(i-1))/(2h) $
e di risolvere il sistema di ODE
$ (d\varphi(t))/dt = aD\varphi $
con un metodo Runge-Kutta su quattro passi il cui array è
$ c = (0,1/2,1/2,1) $
$ A_(ij) = | ( 0 , 0, 0, 0),( 1/2, 0, 0, 0),( 0, 1/2, 0, 0),( 0, 0, 1, 0) | $
$ b = (1/6,1/3,1/3,1/6) $
quindi si ha che
$ \varphi_(n+1) = \varphi_n + dt(1/6K1 + 1/3K2 + 1/3K3 + 1/6K4) $ .
Non ci sono stati problemi nel implementare il metodo RK4. Il problema è che, con quella condizione iniziale, la discretizzazione è instabile.
Ho provato anche a semidiscretizzare la derivata spaziale con una derivata centrata su tre punti, ma ho lo stesso risultato.
L'unica cosa che ho notato è che, se la I.C. è un sinusoide di argomento $ (2\pix) $ allora il tutto funziona e non vi sono problemi.
E' un mio errore nello scrivere la I.C. oppure è lo schema che ha questa "particolarità"?
Ecco il codice
Testo nascosto, fai click qui per vederlo
- Codice:
function K = Kappa1 (phi) %Funzione che traduce il sistema di ODE
global W
phi = phi(:);
K = W*phi;
K = K(:)';
end
function PHI = RK4(t,x,Phi0) %Metodo RK4
dt = t(2) - t(1);
Nt = length(t);
N = length(x);
PHI = nan(Nt,N);
Phi0 = Phi0(:)';
PHI(1,:) = Phi0;
phi = Phi0;
for it=2:Nt
phi1 = phi;
K1 = Kappa1(phi1);
phi2 = phi + (dt/2)*K1;
K2 = Kappa1(phi2);
phi3 = phi + (dt/2)*K2;
K3 = Kappa1(phi3);
phi4 = phi + dt*K3;
K4 = Kappa1(phi4);
phi = phi + dt*(1/6*K1 + 1/3*K2 + 1/3*K3 + 1/6*K4);
PHI(it,:) = phi;
end
end
clc; clear all; close all;
global W
L = 1;
N = 50;
x = linspace (0, L, N+1); x = x(1:end-1); x = x(:);
h = x(2) - x(1);
a = 1;
T = L/a;
Nt = 600;
t = linspace (0, T, Nt);
Phi0 = ((x-L).^2).*cos(2*pi*x/L);
discretizzazione = 1;
switch discretizzazione
case 1
v1 = zeros (1,N);
v2 = zeros (1,N);
v1(1) = 2/3;
v1(2) = 1/6;
v1(end) = 1/6;
A = gallery ('circul', v1);
v2(2) = 1/2;
v2(end) = -1/2;
B = gallery ('circul', (1/h)*v2);
D = (A\B);
W = sparse(N,N);
for i=1:N
W(i,:) = a*D(i,:);
end
case 2
v3 = zeros (1,N);
v3(2) = 1;
v3(end) = -1;
W = gallery ('circul', (a/h)*v3);
end
PHI = RK4(t, x, Phi0);
PHI = [PHI PHI(:,1)];
for i=1:Nt
plot ([x;x(end)+h], PHI(i,:), 'b'); grid on;
drawnow;
end
Grazie.
