Difference between numerical and analytical solution of ode?

3 views (last 30 days)
I G
I G on 21 Apr 2022
Answered: David Goodmanson on 24 Apr 2022
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:
  7 Comments

Sign in to comment.

Answers (1)

David Goodmanson
David Goodmanson on 24 Apr 2022
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
.

Community Treasure Hunt

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

Start Hunting!