# Runge Kutta Fourth order and Interpolation

11 views (last 30 days)
Arunkarthiheyan on 17 Feb 2023
Edited: Arunkarthiheyan on 23 Feb 2023
We are trying to solve a system of coupled differential equations by using RK4 for calculating the values of pressures 'p1' and 'p11'. Using the calculated value of pressure we interpolate the value of Energy density 'ed' in every iteration.
The problem we are facing is that the value of pressure 'p11' is not reducing gradually to reach the surface pressure 'p11(end)'.
Kindly help us rectify the mistake in the code.
clear all
close all
clc
a = importdata('rhoce2.dat') ;
%number density a(:,1) ,energy density a(:,2), parallel pressure a(:,3)
% perpendicular pressure a(:,4)
ed(1) = 5.773955099050900e-06; %intial value for energy density 5899
p11(1) = 4.985073388197410e-10; %intial value for parallel pressure
p1(1) = 5.891730398017590e-10 ; %intial value for perpendicular pressure
r(1) = 1e-6 ; %intial value for perpendicular radius
z(1) = 1e-6 ; %intial value for parallel radius
m1(1)= 4*pi*r(1).^2 *z(1) *ed(1) ; %intial value of perpendicular mass
m11(1) = m1(1) ; %intial value of parallel mass
%p11(end) = 2.5e-19 ; % parallel pressure at surface
%p1(end) = 1.3e-15 ; % perpendicu lar pressure at surface
f1 = @(z,r,ed,p11,m11) -(ed+p11)*z*m11.^3 *((r+1)./(((ed+p11)*(r.^2 +z.^2).^6) -(r/p11))) ./(2*pi*(r.^2 + z.^2)) ; %dp11/dz
f2 = @(r,ed) (4/3)*pi*r.^2 *ed ; %dm11/dz
f3 = @(r,z,ed,p1,m1) -4*pi*(ed+p1)*(r.^2+z.^2).^6 *((ed+p1)-p1*(r-m1)) ./(m1.^3 *r.^2) ; %dp1/dr
f4 = @(r,ed,z) (8/3)*pi*r*ed*z ; %dm1/dr
h = 1e-6 ; % step size for perpendicular radius r
j = 1e-6 ; % step size for parallel radius z
i= zeros(1, 15000)
i=1
%{
for i=1:1000
if p11(i) <=1e-20
break
end
%}
while p11(i)>=1e-20
w1 = h* f1(p11(i),z(i),r(i),m11(i),ed(i)) ; %%%% for z component (first equation) p11
w2 = h* f1(p11(i)+0.5*w1,z(i)+0.5*h,r(i),m11(i)+0.5*w1,ed(i)) ;
w3 = h* f1(p11(i)+0.5*w2,z(i)+0.5*h,r(i),m11(i)+0.5*w2,ed(i)) ;
w4 = h* f1(p11(i)+w3,z(i)+h,r(i),m11(i)+w3,ed(i)) ;
x1 = h* f2(ed(i),r(i)) ; %%%%%%%%% for parallel mass ( 2nd equation) m11
x2 = h* f2(ed(i),r(i)) ;
x3 = h* f2(ed(i),r(i)) ;
x4 = h* f2(ed(i),r(i)) ;
y1 = h* f3(p1(i),z(i),r(i),m1(i),ed(i)) ; %%%% for r component (3rd equation) p1
y2 = h* f3(p1(i)+0.5*y1,z(i),r(i)+0.5*j,m1(i)+0.5*y1,ed(i)) ;
y3 = h* f3(p1(i)+0.5*y2,z(i),r(i)+0.5*j,m1(i)+0.5*y2,ed(i)) ;
y4 = h* f3(p1(i)+y3,z(i),r(i)+j,m1(i)+y3,ed(i)) ;
v1 = h* f4(r(i),ed(i),z(i)) ; %%%%%%%% for perpendicular mass(4th equation) m1
v2 = h* f4(r(i),ed(i),z(i)+0.5*j) ;
v3 = h* f4(r(i),ed(i),z(i)+0.5*j) ;
v4 = h* f4(r(i),ed(i),z(i)+j) ;
p11(i+1) = p11(i) + (1/6)* (w1+2*w2+2*w3+w4);
m11(i+1) = m11(i) + (1/6)* (x1+2*x2+2*x3+x4);
p1(i+1) = p1(i) + (1/6)* (y1+2*y2+2*y3+y4);
m1(i+1) = m1(i) + (1/6)* (v1+2*v2+2*v3+v4);
r(i+1) = r(i) + h ;
z(i+1) = z(i) + j ;
%ed(2)= interp1( a(:,2), a(:,4) , p1(2), 'linear');
%ed(i+1)= interp1(a(:,2), p1(i+1));
% ed(i+1)= griddedInterpolant((a(:,2)), (a(:,4)), linear, p1(i+1) );
%F= griddedInterpolant(unique (a(:,4)) , unique ( a(:,2)) );
%ed(i+1) = F(p1(i+1)) ;
%%%interpolation
en = unique ( a(:, 2) ) ; % energy density
pnp = unique ( a(:, 4) ) ; % perpendicular pressure
%p1(2)= 5.790305785983432e-10;
x11= 0; x22= 0; y11= 0; y22= 0;
for k= 2 : 14601
if (p1(i+1) < pnp(k))
x11= en(k-1) ;
x22= en(k);
y11= pnp(k-1);
y22= pnp(k);
break
end
end
ed(i+1) = x11 + (x22-x11) * ((p1(i+1)-y11)/(y22-y11)) ;
i = i+1
end
Torsten on 17 Feb 2023
You have differential equations in r and differential equations in z. You cannot mix them in a "combined" Runge-Kutta-method as you obviously try to do in your code.
Kishan Bajpai on 18 Feb 2023
I try in that way , Thanks

Alan Stevens on 17 Feb 2023
Edited: Alan Stevens on 17 Feb 2023
Shouldn't
f1 = @(z,r,ed,p11,m11) -(ed+p11)*z*m11.^3 *((r+1)./(((ed+p11)*(r.^2 +z.^2).^6) -(r/p11))) ./(2*pi*(r.^2 + z.^2)) ;
be
f1 = @(z,r,ed,p11,m11) -(ed+p11)*z*m11.^3 *((r.^2+r)./(((ed+p11)*(r.^2 +z.^2).^6) -(r/p11))) ./(2*pi*(r.^2 + z.^2)) ;
Can't check anything much as you don't supply file rhoce2.dat.
Arunkarthiheyan on 22 Feb 2023
Edited: Arunkarthiheyan on 23 Feb 2023
I checked the code after changing the 'f1' which you have indicated, but still the error persists.
rhoce2 file has been attached herewith.
Kindly help us rectify the problem. Thanks @Alan Stevens