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 Center および File ExchangeOrdinary Differential Equations についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by