Esercizio BDF2
Ciao, ho questo problema di diffuzione-trasporto
\(\displaystyle du/dt=1/100 d^2u/dx^2 +6du/dx\)
\(\displaystyle u(0,x)=5x^2(x-4/5)\)
\(\displaystyle du/dx(t,0)=0\)
\(\displaystyle u(t,1)=1\)
Devo risolverlo usando differenze finite centrate nello spazio e mostrare il corretto ordine di convergena del metodo BDF di ordine 2 al tempo t*=0.1
Io ho scritto con matlab questo codice
Compilandolo mi da errore nel calcolo di y(2), ma non so come sistemare il codice.
Qualcuno ha una possibile soluzione?
Grazie
\(\displaystyle du/dt=1/100 d^2u/dx^2 +6du/dx\)
\(\displaystyle u(0,x)=5x^2(x-4/5)\)
\(\displaystyle du/dx(t,0)=0\)
\(\displaystyle u(t,1)=1\)
Devo risolverlo usando differenze finite centrate nello spazio e mostrare il corretto ordine di convergena del metodo BDF di ordine 2 al tempo t*=0.1
Io ho scritto con matlab questo codice
clear all close all clc a = 0; b = 1; c = 6; d = 1/100; % differenze finite nello spazio m = 51; x = linspace(a,b,m)'; h = (b-a)/(m-1); A = d * toeplitz(sparse([1,1],[1,2],[-2,1]/h^2,1,m)); A(1,1:2) = d * [-2/h^2,2/h^2]; B = c * toeplitz(sparse (1,2,-1/(2*h),1,m), sparse(1,2,1/(2*h),1,m)); B(1,1:2) = [0,0]; C = d*A+c*B; C(m,m-1:m) = [0,0]; Peclet = abs(c)*h/(2*d) % uso un metodo di ordine 2 nel tempo tstar = 0.1; ts = 700; k = tstar/ts; y0 = @(x) 5*x.^2.*(x-4/5); yn = y0; Yes=y0(x); tsrange = [10 20 50 100 200]; counter = 0; for ts = tsrange counter = counter+1; k = tstar/ts; y(1) = y0; y(2) = (1-(k/2)*C) \ (1+(k/2)*C)*y(1); % trapezi for n = 1:ts-1 y(n+2)= (1-(2/3)*k*C)\(-(1/3)*y(n)+(4/3)*y(n+1)); end errore(counter) = norm(Yes-y,inf); L(counter) = length(x); end loglog(L,errore,'-*',L,L.^(-1),'g',L,L.^(-2),'r',L,L.^(-3),'b',L,L.^(-4),'y') legend('errore del metodo','retta pendenza -1','retta pendenza -2','retta pendenza -3','retta pendenza -4'); title('Ordine di convergenza del metodo'); yrif = yn;
Compilandolo mi da errore nel calcolo di y(2), ma non so come sistemare il codice.
Qualcuno ha una possibile soluzione?
Grazie
Risposte
Se per calcolare $y(2)$ devi fare una divisione allora devi usare /, non \.