Difference between numerical and analytical solution of ode?

29 views (last 30 days)
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
.

Categories

Find more on Numerical Integration and 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!

Translated by