Usually, the numerical solvers in Matlab work well with the standard settings. Sometimes they do not, and it is not always obvious they have not worked! Part of using a tool like Matlab is checking how well your solution really worked. We use an example of integrating an ODE that defines the van der Waal equation of an ideal gas here.
clc; clear all; close all
we plot the analytical solution to the van der waal equation in reduced form here.
Tr = 0.9; Vr = linspace(0.34,4,1000); % analytical equation for Pr Prfh = @(Vr) 8/3*Tr./(Vr - 1/3) - 3./(Vr.^2); Pr = Prfh(Vr); % evaluated on our reduced volume vector. % Plot the EOS plot(Vr,Pr) ylim([0 2]) xlabel('V_R') ylabel('P_R')
we want an equation for dPdV, which we will integrate we use symbolic math to do the derivative for us.
syms Vrs Trs Prs = 8/3*Tr./(Vrs - 1/3) - 3./(Vrs.^2) ; diff(Prs,Vrs)
ans = 6/Vrs^3 - 12/(5*(Vrs - 1/3)^2)
Vspan = [0.34 4]; Po = Prfh(Vspan(1)); [V,P] = ode45(@myode,Vspan,Po); hold all plot(V,P)
That is strange, the numerical solution does not agree with the analytical solution. Let's try a few more solvers:
[V,P] = ode15s(@myode,Vspan,Po); plot(V,P) [V,P] = ode113(@myode,Vspan,Po); plot(V,P) legend 'analytical' 'ode45' 'ode15s' 'ode113'
None of the solvers give good solutions! Let's take a look at the problem. First, look at the derivative values
figure; hold all plot(Vr,myode(Vr,P)) % analytical derivative plot(V,cmu.der.derc(V,P)) % numerical derivative xlabel('V_R') ylabel('dPdV') % these look the same, but let's note the scale is 10^4! Lets look at % relative errors:
dPdV = cmu.der.derc(V,P); error = (dPdV - myode(V,P))./dPdV; figure plot(V,error) xlabel('V_R') ylabel('relative error') ylim([-0.1 0.1]) % you can see here the errors in the derivative of the numerical solution % are somewhat large, sometimes up to 5%. That suggests the tolerances are % not set appropriately for the ODE solver, or that what is accurate for % 10^4, is not also accurate when dPdV is very small.
We tighten the tolerances from 1e-3 (the default) to 1e-9
options = odeset('AbsTol',1e-9,'RelTol',1e-9); [V,P] = ode45(@myode,[0.34,4], Prfh(0.34),options); figure hold all plot(Vr,Pr) plot(V,P,'r--') ylim([0 2]) xlabel('V_R') ylabel('P_R')
you can see much closer agreement with the analytical solution now.
Now, the errors are much smaller, with a few exceptions that don't have a big impact on the solution.
dPdV = cmu.der.derc(V,P); error = (dPdV - myode(V,P))./dPdV; figure plot(V,error) xlabel('V_R') ylabel('relative error') ylim([-0.1 0.1])
The problem here was the derivative value varied by four orders of magnitude over the integration range, so the default tolerances were insufficient to accurately estimate the numerical derivatives over that range. Tightening the tolerances helped resolve that problem. Another approach might be to split the integration up into different regions. For instance, if instead of starting at Vr = 0.34, which is very close to a sigularity in the van der waal equation at Vr = 1/3, if you start at Vr = 0.5, the solution integrates just fine with the standard tolerances.
ans = done
function dPrdVr = myode(Vr,Pr) dPrdVr = 6./Vr.^3 - 12./(5*(Vr - 1/3).^2); % categories: odes