Difference between numerical and analytical solution of ode?
7 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/972555/image.jpeg)
7 commentaires
Réponses (1)
David Goodmanson
le 24 Avr 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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/976100/image.png)
% 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 commentaires
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!