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

1 view (last 30 days)
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

Answers (0)

Community Treasure Hunt

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

Start Hunting!