matlab - كيفية الحصول على دقة فورتران في ماتلاب



fortran double-precision (1)

كنت لا تزال لا تستخدم بالتالي دقة مزدوجة في جميع أنحاء التعليمات البرمجية الخاصة بك، على سبيل المثال:

beta     = 12.D0*0.0001/(1.D0*( (1.0 - 0.1)**4 ))

و أكثر من ذلك بكثير. إذا كنت إجبار المترجم على استخدام الدقة المزدوجة كما الافتراضي ل gfortrangfortran الخيار تجميع هو -fdefault-real-8 )، والنتيجة من التعليمات البرمجية الخاصة بك:

0،00000000000000000000000000000000000

لذلك تحتاج إلى إصلاح التعليمات البرمجية الخاصة بك. وعلى سبيل المثال، ينبغي أن يقرأ السطر المذكور:

beta     = 12.D0*0.0001D0/(1.D0*( (1.0D0 - 0.1D0)**4 ))

[على الرغم من أنني يحتقر تدوين D0 ، ولكن هذا هو قصة مختلفة]

لدي قطعة من التعليمات البرمجية مكتوبة في فورتران وفي ماتلاب. أنها تفعل بالضبط نفس الحساب، وهي

  1. بناء حقل tanh والعثور على لابلاسيان لها
  2. مضاعفة بعض المصطلحات معا

نتيجة هذا الضرب ينتج مصفوفة، الذي (4،4) ث و (6،6) ط طرح.

  • في فورتران الفرق بينهما هو ~ 1e-20
  • في ماتلاب الفرق بينهما هو صفر تقريبا.

هذه المسألة حرجة جدا، وأنا اختبر إذا كان هذا الرقم أقل من الصفر. سؤال : هل هناك طريقة لأداء الحسابات بحيث أحصل على نفس الدقة في ماتلاب كما في فورتران؟

أدرج الرموز أدناه:

MATLAB

clear all

weights = [4./9, 1./9,1./9,1./9,1./9, 1./36,1./36,1./36,1./36];
dir_x   = [  0,   1,  0, -1,  0,    1,  -1,  -1,   1];
dir_y   = [  0,   0,  1,  0, -1,    1,   1,  -1,  -1];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CONSTANTS
length_y = 11; length_x = length_y;
y_center = 5; x_center  = y_center;


densityHigh = 1.0;
densityLow  = 0.1;
radius  = 3.0;
c_width = 1.0;

average_density = 0.5*(densityHigh+densityLow);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





for x=1:length_x
    for y=1:length_y
        for i=1:9
            fIn(i, x, y) = weights(i)*densityHigh;
            test_radius = sqrt((x-x_center)*(x-x_center) + (y-y_center)*(y-y_center));
            if(test_radius <= (radius+c_width))
                fIn(i, x, y) = weights(i)*( average_density - 0.5*(densityHigh-densityLow)*tanh(2.0*(radius-sqrt((x-x_center)*(x-x_center) + (y-y_center)*(y-y_center))/c_width)) );
            end
        end
    end
end 

ref_density_2d = ones(length_x)*average_density;
for i=1:length_x
    ref_density(:,:,i) = abs(ref_density_2d(:, i)');
end


          rho = sum(fIn);
laplacian_rho = (+1.0*(circshift(rho(1,:,:), [0, -1, -1]) + circshift(rho(1,:,:), [0, +1, -1]) + circshift(rho(1,:,:), [0, -1, +1]) + circshift(rho(1,:,:), [0, +1, +1])) + ...
                 +4.0*(circshift(rho(1,:,:), [0, -1, +0]) + circshift(rho(1,:,:), [0, +1, +0]) + circshift(rho(1,:,:), [0, +0, -1]) + circshift(rho(1,:,:), [0, +0, +1])) + ...
                -20.0*rho(1,:,:));

psi   = 4.0*0.001828989483310*(rho-densityLow).*(rho-densityHigh).*(rho-ref_density) - laplacian_rho*(1.851851851851852e-04)/6.0;

psi(1,4,4)-psi(1,6,6) 

فورتران

 PROGRAM main

 IMPLICIT NONE

 INTEGER, PARAMETER :: DBL = KIND(1.D0)
 REAL(KIND = DBL), DIMENSION(1:11,1:11) :: psi, rho
 INTEGER :: i, j, m, ie, iw, jn, js
 REAL(KIND = DBL) :: R, rhon, lapRho

 INTEGER, DIMENSION(1:11,1:11,1:4) :: ni

 REAL(KIND = DBL) :: kappa, kappa_6, kappa_12, kappaEf, beta, beta4



 beta     = 12.D0*0.0001/(1.D0*( (1.0 - 0.1)**4 ))
 kappa    = 1.5D0*0.0001*1.D0/( (1.0 - 0.1)**2 ) 


!-------- Define near neighbours and initialize the density rho ----------------
 DO j = 1, 11
   DO i = 1, 11

! Initialize density
      rho(i,j) = 1.D0
        R =  DSQRT( ( DBLE(i)-5.0 )**2 + ( DBLE(j)-5.0 )**2 )
        IF (R <= (DBLE(3.0) + 1.D0)) THEN
          rho(i,j) = 0.55D0 - 0.5*0.9*TANH(2.D0*(DBLE(3.0) - R)/1.D0)
        END IF

 !Generate neighbors array
      ni(i,j,1) = i + 1
      ni(i,j,2) = j + 1
      ni(i,j,3) = i - 1
      ni(i,j,4) = j - 1
   END DO
 END DO


! Fix neighbours at edges
 ni(1,:,3) = 11
 ni(11,:,1) = 1
 ni(:,1,4) = 11
 ni(:,11,2) = 1



!--------- Differential terms for the stress form of the interfacial force -----
 DO j = 1, 11
   DO i = 1, 11

! Identify neighbors
     ie = ni(i,j,1)
     jn = ni(i,j,2)
     iw = ni(i,j,3)
     js = ni(i,j,4)

! Laplacian of the density rho
     lapRho = 4.D0*( rho(ie,j ) + rho(iw,j ) + rho(i ,jn) + rho(i ,js) )       &
             + rho(ie,jn) + rho(ie,js) + rho(iw,jn) + rho(iw,js) - 20.D0*rho(i,j)

! Define the chemical potential Psi
     psi(i,j) = 4.D0*beta*( rho(i,j) - 0.55 )*( rho(i,j) - 0.1 )*( rho(i,j) - 1.0 ) &
              - kappa*lapRho/6.D0
   END DO
 END DO


write(*,*) psi(6,6)-psi(4,4)


 END PROGRAM




double-precision