Comparing an analytic solution to a numerical solution (it does not work!)

조회 수: 1 (최근 30일)
Douglas Alves
Douglas Alves 2014년 5월 14일
There`s an equation called filament solution by Ostriker. I need to compare his solution to the solution which comes from a system of diferential equation I`m supposed to find a ratio of almost zero. (1e-17) but both curves seems to be the same but they don`t completely match. Am I doing something wrong?
function poisson_values
par=setup;
inits = [1 0 0 ];
Rrange = [1e-03 1e+02];
options = odeset('RelTol',1e-4,'AbsTol',1e-5);
[r, y] = ode23t(@p,Rrange,inits,options,par);
figure(1);
hold on
loglog(r,y(:,1),'g'); grid; xlabel('log(r)') ; ylabel('log(rho)');
Solut(r);
ratio = ( output - y(:,1));
%x_axis = 1:149;
%ratio1= abs(ratio);
figure(2);
plot(abs(ratio),'r'); ylabel('ratio') ; grid ;
legend('ratio');
function Ost_solution = Solut(r,par)
par=setup;
r0 = (par.sigma)./sqrt(4*pi*par.G*par.rho_c);
Ost_solution = (par.rho_c)./(1+(r.^2)./8*(r0.^2)).^2 ;
output = Ost_solution ;
loglog(r,Ost_solution,'b') ; xlabel('log(r)') ; ylabel('log(rho)');
legend('ode45 solver', 'analytic solution');
hold off
end
end
function poisson = p(r,y,par)
rho = y(1);
m = y(2);
g=y(3);
%I = y(4);
poisson = [ -rho.*g;2*pi.*r*rho;rho - g./r];
end
function par=setup
par.rho_c = 1.0 ;
par.G = 1.0;
par.sigma = 1.0 ;
end

답변 (0개)

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by