Difference between numerical and analytical solution of ode?

I am solving two differential equations numericaly and analyticaly. But when I compare results, they are not the same, I cannot conclude why?

These are equations:

x0' = - 32 .* beta ./ (x0 .* R .^ 4)
x1' = (- 8 .* x0' ./ R - x0' .* x1) ./ x0

Also, I have next conditions:

z=1: x0=1, x1=0

And

x0' = d (x0) / dz, R = R(z) = Ri - z .* (Ri - 1)

When I solve it analyticaly, I got:

x0 = (1 + 64 .* beta .* (1 - 1 ./ R .^ 3) ./  (3 .* (Ri - 1))) .^ 0.5
x1 = 8 .* (1 .* x0 - 1) ./ R

I am solving it numericaly in this way:

function [f, R] = fun_p(z, x, beta, ri)
R = ri - z .* (ri - 1);
f = zeros(2, size(x,2));
f(1,:) = - 32 .* beta ./ (R .^ 4 .* x(1,:));
f(2,:) = ( - 8 .* f(1,:) ./ R - f(1,:) .* x(2,:) ) ./ x(1,:);

I am calling fun_p from this file:

clear all;
clc;
options = odeset('RelTol',1.e-6, 'AbsTol',1.e-6);
Ree = 0.1;
Kne = 0.1;
eps = 0.01;
beta = 0.06;
z = linspace(1, 0, 1001);
ri = 0.7;
R = ri - z .* (ri - 1);
[~, pv] = ode45(@(z, x)fun_p(z, x, beta, ri), z, [1; 0], options);
x0 = pv(:, 1);
x1 = pv(:, 2);
x00 = ( 1 + 64 .* beta .* ( 1 - 1 ./ R .^ 3) ./ ( 3 .* (ri - 1))) .^ 0.5;
x11 = 8 .* ( 1 ./ x00 - 1 ) ./ R;

figure;
plot(z,x00);
hold on;
plot(z,x0, 'x');
hold on;

plot(z,x11);
hold on;
plot(z,x1, 'x');
hold on;
legend({'analytical', 'numerical', 'analytical', 'numerical'}, 'FontSize', 16, 'Interpreter', 'LaTeX');

When I compare results I get som difference for x1, I cannot conclude why:

[Torsten]

Maybe it's trivial, but even with symbolic computation, I can't get

x1' - (- 8 .* x0' ./ R - x0' .* x1) ./ x0 == 0

for your analytical "solutions" x0 and x1.

[I G]

@Torsten I wrote it in this form:

x1' = (- 8 .* x0' ./ R - x0' .* x1) ./ x0 /multiply with x0

x1' * x0 = - 8 .* x0' ./ R - x0' .* x1

x1' * x0 + x0' .* x1 = - 8 .* x0' ./ R

(x1 * x0)' = - 8 .* x0' ./ R

So, I am solving it symbolic as:

d( x0 *x1 ) / dz = ( - 8 * dx0 / dz) / R /multiply with dz

d(x0 * x1) = - 8 * dx0 / R / integrate

x0 * x1 = - 8 * x0 / R + C

z = 1: x0 = 1, x1 = 0

0 = -8 / R + C, C = 8 / R

x0 * x1 = - 8 * x0 / R + 8 / R

x1 = 8 * (1 / x0 - 1) / R

I am just confused with this, R is dependent on z: R = ri - z .* (ri - 1). Should I take this dependency into account somewhere in this integration, or in whole equation multiplication with dz? And how to take it into account and integrate in that case?

[I G]

I took your advice, and this is now my solution procedure:

d( x0 * x1 ) / dz = ( - 8 * dx0 / dz) / R

d( x0 * x1 ) = - 8 * ( d( x0 / R ) / dz) dz /*how did you pulled in R into d,so dx0 -> d(x0/R)

d( x0 * x1 ) = - 8 * (dx0/dz / R - x0/R^2 dR/dz) dz /integrate

x0 * x1 = - 8 * (x0 / R - x0/R^2 ( - ( ri-1 )) * z) + C

z = 1: x0 = 1, x1 = 0

0 = - 8 * (1 / R + 1 / R^2 (ri-1)) + C

C = 8 * (1 / R + 1 / R^2 (ri-1))

x0 * x1 = - 8 * (x0 / R + x0/R^2 ( ri-1 ) * z) + 8 * (1 / R + 1 / R^2 (ri-1))

x1 = 8 * (- 1 + 1 / x0) / R + 8 * (ri - 1) * (1 / x0 - z) / R^2

When I plot x1, and compare with numerical solution with ode45, it is still not the same, I don`t know why? x0 is the same as the numerical solution.

[I G]

@Bruno Luong How did you pulled in R into d, in previous answer:

(( - 8 * dx0 / dz) / R) dz -> -8 * d(x0/R)/dz dz

Is that possible, when in original equation it is only (dx0/dz)/R, how did you involve R into this derivative d(x0/R)/dz?

And what is correct solution for x1 according to you?

[I G]

Well I don't ask for solutions, I ask for instructions. I wouldn't be here if it was completely clear to me.

[I G]

No, I don`t understant how did you involved R into differential (into f(z) in your example), when at the start there is (1/R(z))*(dx0/dz)? How is that: (1/R(z))*(dx0/dz) = d(x0/R(z))/dz

Also, I don`t understand, you said "No wrong integral fix.", but there is still some difference between analytical and numerical solution, where I should look further? Or on which other way I can try to solve x1?

[Torsten]

I doubt there is an analytical solution for x1.

Answer 1

[David Goodmanson]

Hi IG,

Here is an almost-analytical solution for x0 and x1 (which are called x and y respectively). It's analytic for x, uses numerical integration for y (not the same integration that is effectively used by ode45) and agrees with the ode45 results.

% x' = - 32*beta/(x*R^4)
% y' = (- 8 .* x' ./ R - x' .* y) ./ x   % (1)
% @ z=1: x0=1, x1=0

beta = .06;
Ri = .7;

[z u] = ode45(@(z,u) fun(z,u,beta,Ri), [1 0],[1 0]);
x = u(:,1);
y = u(:,2);

% analytic and numerical integration solution
% (x*y)' = 256*beta/(x*R^5);     % equivalent to (1) 
za = linspace(0,1,1000);
R = Ri - za*(Ri - 1);
xa = sqrt(-64*beta./(3*(Ri-1)*R.^3) + 64*beta/(3*(Ri-1)) +1 );
I = cumtrapz(za,256*beta./(xa.*R.^5));
ya = (1./xa).*(I-I(end));

fig(2)
plot(z,x,'o',za,xa,z,y,'o',za,ya)
legend('x ode45','x analytic','y ode45','y by integration','location', 'southeast')

function dxydz = fun(z,u,beta,Ri)
R = Ri - z*(Ri - 1);
x = u(1);
y = u(2);
dxdz = -32*beta/(x*R^4);
dydz = -(dxdz/x)*(8/R +y);      % (1)
dxydz = [dxdz; dydz];
end  

.