3 added original equations
source | link

The equations I want to solve are as below. I want to solve v(z),qr(z) and qi(z) and plot v(z0)&b(z0)&qi(z0) v.s.z0, where b=Sqrt(1/qr)/k0, z=z0*k0. The initial conditions are, v(0)=1,b0=100,qi(0)=0, N=100.

enter image description here

enter image description here

enter image description here The plots from mma(left) and matlab(right) are: enter image description here

enter image description here

enter image description here

The equations I want to solve are as below. I want to solve v(z),qr(z) and qi(z) and plot v(z0)&b(z0)&qi(z0) v.s.z0, where b=Sqrt(1/qr)/k0, z=z0*k0. The initial conditions are, v(0)=1,b0=100,qi(0)=0, N=100.

enter image description here

The plots from mma(left) and matlab(right) are: enter image description here

2 added Rungekutta method in NDSolve and graphs in mma and matlab
source | link

I am solving nonlinear ODE + integration groups using NDSolve of mma and ode45 of matlab. At last, I found that they give different results, even after I added RungeKutta method, interpolationorder, changed the value of precisiongoal, and accuracygoal, methods,interpretionorder of NDSolve. Has anyone met similar problemsDo the different answer origin from numerical error, or the wrong coding? 

MMA codes and results (solve v(z),qr(z),qi(z); b=1/qr^2; plot v(z0), b(z0), qi(z0) v.s. z0=z/k0)

enter image description here

Matlab codes and results, y(1)=v, y(2)=qr, y(3)=qi:

enter image description here

I am solving nonlinear ODE + integration groups using NDSolve of mma and ode45 of matlab. At last, I found that they give different results, even after I changed the value of precisiongoal, accuracygoal, methods,interpretionorder of NDSolve. Has anyone met similar problems?

MMA codes (solve v(z),qr(z),qi(z); b=1/qr^2; plot v(z0), b(z0), qi(z0) v.s. z0=z/k0)

Matlab codes, y(1)=v, y(2)=qr, y(3)=qi:

I am solving nonlinear ODE + integration groups using NDSolve of mma and ode45 of matlab. At last, I found that they give different results, even after I added RungeKutta method, interpolationorder, changed the value of precisiongoal, and accuracygoal of NDSolve. Do the different answer origin from numerical error, or the wrong coding? 

MMA codes and results (solve v(z),qr(z),qi(z); b=1/qr^2; plot v(z0), b(z0), qi(z0) v.s. z0=z/k0)

enter image description here

Matlab codes and results, y(1)=v, y(2)=qr, y(3)=qi:

enter image description here

1
source | link

Different results from ode45 and NDSolve?

I am solving nonlinear ODE + integration groups using NDSolve of mma and ode45 of matlab. At last, I found that they give different results, even after I changed the value of precisiongoal, accuracygoal, methods,interpretionorder of NDSolve. Has anyone met similar problems?

MMA codes (solve v(z),qr(z),qi(z); b=1/qr^2; plot v(z0), b(z0), qi(z0) v.s. z0=z/k0)

Clear["`*"]
np = 8.;
\[Beta] = 9.5*10^-7;
k0 = 2 Pi;
I0I14 = 1.;
b0 = 100.;
z0max = 100000.;
zmax = z0max*k0;
qr0 = 1/(k0*b0)^2;

ek = 6*10^-5*I0I14* v^2*Exp[-2*qr*rho^2];
ep = 1.8*10^-4*1*I0I14^8*(v^2*Exp[-2*qr*rho^2])^8;
er = 1. + ek + ep;
ei = 9.5*10^-7*I0I14^7* (v^2*Exp[-2*qr*rho^2])^7;

{coreA[v_, qr_, qi_, rho_], coreB[v_, qr_, qi_, rho_]} = 
Exp[-qr rho^2] rho RotationMatrix[-qi rho^2].{ei, 1 - er};
int[core_, v_, qr_, qi_?NumericQ] := 
v NIntegrate[core[v, qr, qi, rho], {rho, 0, Infinity}, MaxRecursion -> 40]

intD = \[Beta]/np*I0I14^7;
eqn = With[{v = v@z, qr = qr@z, qi = qi@z}, 
With[{intA = int[coreA, v, qr, qi], intB = int[coreB, v, qr, qi]},
{D[v, z] == -intA qr - intD v^(2*np - 1),
 D[qr, z] == -2 intA qr^2/v - intD qr v^(2*np - 2),
 D[qi, z] == -3 intA qi qr/v + intB qr^2/v - intD qi v^(2*np - 2)
 }]];

bc = {v[0] == 1., qr[0] == qr0, qi[0] == 0.};

sol = NDSolve[{eqn, bc}, {v, qr, qi}, {z, 0., zmax}, InterpolationOrder -> All] 

LogLinearPlot[v[z*k0] /. sol[[2]] // Evaluate, {z, z0min, z0max}]
LogLinearPlot[{Sqrt[1/qr[z*k0]]/k0} /. sol[[2]] // Evaluate, {z, z0min, z0max}]
LogLinearPlot[{qi[z*k0]} /. sol[[2]] // Evaluate, {z, z0min, z0max}]

Matlab codes, y(1)=v, y(2)=qr, y(3)=qi:

clear all

k0 = 2 * pi;    nu_0 = 1;                     
b0 = 100;                        
b_0 = b0 * k0;                  
qr_0 = 1/b_0^2;                    % initial qr
qi_0 = 0;                         

options = odeset('RelTol',1e-3,'AbsTol',[1e-4 1e-4 1e-4]);
ts0 = 0; 
tf0 = 1E5;                         
ts = round(0 * k0); 
tf = round(tf0 * k0);
nsteps = 1E5;
step = (tf - ts) / nsteps;
tspan = linspace(ts, tf, nsteps);

T = zeros(nsteps,1);
F = zeros(nsteps,3);
warning('OFF','MATLAB:integral:MinStepSize');
[T,F] = ode45(@test_fun,tspan,y0,options);

%------------Plot--------------%

z0 = T / k0;
nu = F(:,1);
qr = F(:, 2);
qi = F(:, 3);
b0 = 1./(F(:,2)).^0.5 / k0;

plot(z0, nu);      
set(gca,'XScale','log')
plot (z0, b0);     
set(gca,'XScale','log')
plot (z0, qi);    
set(gca,'XScale','log')

-----------function------------
function dy=test_fun(x,y)
I14 = 1;
beta = 9.5 * 10^(-7);
np = 8;

ee = @(r)(I14 * y(1)^2 * exp(-2 * y(2) * r.^2));
ek = @(r) (6 * 10^(-5) * ee(r));
ep = @(r) (1.8 * 10^(-4) * 1 * ee(r).^8);
er = @(r) (1 + ek(r) + ep(r));
ei = @(r) (9.5 * 10^(-7) * ee(r).^7);

% A, B, D
A_func = @(r)(exp(-y(2) * r.^2).* ((1 - er(r)).* sin(y(1) * r.^2 * pi/180) + ei(r).* cos(y(3)* r.^2 * pi/180)).* r);
A = y(1) * integral(A_func, 0, Inf);

B_func = @(r)(exp(-y(2) * r.^2).* ((1 - er(r)).* cos(y(1) * r.^2 * pi/180) - ei(r).* sin(y(3) * r.^2 * pi/180)).* r);
B = y(1) *integral(B_func, 0, Inf);

D = beta * I14^7 / np;


dy = zeros(3,1);
dy(1) = -A * y(2) - D * y(1)^(2 * np - 1);
dy(2) = -2 * A * y(2)^2 / y(1) - D * y(2) * y(1)^(2 * np -2);
dy(3) = -3 * A * y(2) * y(3) / y(1) + B * y(2)^2 / y(1) - D * y(3) * y(1)^(2 * np - 2);