Comparing an analytic solution to a numerical solution (it does not work!)
1 view (last 30 days)
Show older comments
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 Comments
Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!