Difference between numerical and analytical solution of ode?
5 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
David Goodmanson
am 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
.
0 Kommentare
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!