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

1 次查看（过去 30 天）

Douglas Alves 2014-5-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/129657-comparing-an-analytic-solution-to-a-numerical-solution-it-does-not-work

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