05/07/2019, 16:03
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
05/07/2019, 18:56
D1 = toeplitz(sparse(1,2,-1/(2*h),1,m),sparse(1,2,1/(2*h),1,m));
05/07/2019, 22:13
feddy ha scritto:Ho letto in fretta... come hai imposto le condizioni al bordo periodiche?
Ad ogni modo, questo è il metodo delle linee... discretizzi gli operatori spaziali (in questo caso con un'accuratzza al secondo ordine) e poi risolvi un sistema di ODEs. Nota che visto che la PDE di partenza è lineare, anche la ODE sarà lineare. E' un semplice problema del trasporto (scalare). La matrice che discretizza la derivata prima tra l'altro si può costruire in un colpo solo come ho già fatto in altre discussioni, così ad esempio:
- Codice:
D1 = toeplitz(sparse(1,2,-1/(2*h),1,m),sparse(1,2,1/(2*h),1,m));
dove $h$ è l'ampiezza dell'intervallo, m il numero di nodi.
Non ho tempo di guardare la tua implementazione di RK4... la suppongo corretta comunque
Non è che hai preso un passo troppo largo?
v1 = zeros (1,N);
v2 = zeros (1,N);
v1(1) = 2/3;
v1(2) = 1/6;
v1(end) = 1/6;
v2(2) = 1/2;
v2(end) = -1/2;
A = gallery ('circul', v1);
B = gallery ('circul', (1/h)*v2);
D = (A\B);
W = sparse(N,N);
for i=1:N
W(i,:) = a*D(i,:);
end
06/07/2019, 09:56
06/07/2019, 10:15
feddy ha scritto:Dove vedi l'operatore di derivata seconda?
11/07/2019, 20:10
11/07/2019, 20:28
Skuola.net News è una testata giornalistica iscritta al Registro degli Operatori della Comunicazione.
Registrazione: n° 20792 del 23/12/2010.
©2000—
Skuola Network s.r.l. Tutti i diritti riservati. — P.I. 10404470014.
Powered by phpBB © phpBB Group - Privacy policy - Cookie privacy
phpBB Mobile / SEO by Artodia.