Runge Kutta Fourth order and Interpolation

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]

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]

I try in that way , Thanks

Answer 1

[Alan Stevens]

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]

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