How to convert fortran code to matlab?

3 次查看（过去 30 天）
Alfandias Saurena
Alfandias Saurena 2022-3-10
评论： Ben Barrowes 2023-7-25
PROGRAM IMPLISIT
C
C *** PROGRAM PENYELESAIAN PERSAMAAN HYPERBOLIK DENGAN METODA
IMPLICIT
C
PARAMETER( IMX=1000 )
DIMENSION U(IMX),UL(IMX),UB(IMX),XP(IMX),UP(IMX,100),TP(IMX),
: A(IMX),B(IMX),C(IMX),D(IMX)
C
C-------------------------------------------------------------------
-----------
C
C UL(I) = HARGA LAMA VARIABLE PADA TITIK I
C UB(I) = HARGA BARU VARIABLE PADA TITIK I
V-12
C UP(I,K) = BEBERAPA HARGA VARIABLE PADA TITIK I DISIMPAN UNTUK PLOTTING
C
C------------------------------------------------------------------------------
C
OPEN(2,FILE='IM1.CON')
OPEN(3,FILE='IM1.OUT')
C
C *** INPUT DATA DAN PARAMETER
C
IM = 801
XL = 0
XR = 20
DX = (XR - XL )/(IM-1)
DT = 0.5
TMAX = 101
CR = 0.1
THETA2 = 0.5
THETA1 = 1.0 - THETA2
IPRINT = 20
CN = CR*DT/DX
C
C *** KOORDINAT NODE
C
DO I=1,IM
XP(I) = XL + (I-1)*DX
ENDDO
C
C *** HARGA AWAL (INITIAL VALUES)
C
DO I=1,IM
IF( XP(I).GT.2.5.AND.XP(I).LE.3.0 )THEN
UB(I) = ( XP(I) - 2.5 )
ELSEIF( XP(I).GT.3.0.AND.XP(I).LE.3.5 )THEN
UB(I) = ( 3.5 - XP(I) )
ELSE
UB(I) = 0.0
ENDIF
ENDDO
C
C *** TETAPKAN HARGA UB(I) PADA SYARAT BATAS EXTERNAL
C
UB(1) = 0
UB(IM) = 0
C
C *** LOOP ITERASI FTCS
C
ITER = 0
IP = IPRINT
TIME = 0
K = 0
100 CONTINUE
V-13
C
C ***** SIMPAN HARGA UB(I) PADA TIME STEP N
C
DO I=1,IM
UL(I) = UB(I)
ENDDO
C
C **** SIMPAN HASIL UNTUK PLOTTING
C
IF( IP.EQ.IPRINT )THEN
IP = 0
K = K + 1
DO I=1,IM
UP(I,K) = UB(I)
TP(K) = TIME
ENDDO
WRITE(*,'(I10,3F10.3)') ITER,TIME,CN
WRITE(2,'(I10,100F10.3)') ITER,TIME,(UB(I),I=1,IM)
ENDIF
ITER = ITER + 1
TIME = TIME + DT
IP = IP + 1
C
C ***** PEMBENTUKAN MATRIK TRIDIAGONAL UNTUK MEMPERBAHARUI HARGA UB(I) PADA TIME STEP N+1
C
DO I=2,IM-1
II = I - 1
IF( I.EQ.2 )THEN
A(II) = 0
B(II) = 1.0
C(II) = THETA1*CN/2
D(II) = UL(I) + THETA2*CN/2*( UL(I+1) - UL(I-1) ) +
: THETA1*CN/2*UL(I-1)
ELSEIF( I.EQ.IM-1 )THEN
A(II) = -THETA1*CN/2
B(II) = 1.0
C(II) = 0
D(II) = UL(I) - THETA2*CN/2*( UL(I+1) - UL(I-1) ) -
: THETA1*CN/2*UL(I+1)
ELSE
A(II) = -THETA1*CN/2
B(II) = 1.0
C(II) = THETA1*CN/2
D(II) = UL(I) - THETA2*CN/2*( UL(I+1) - UL(I-1) )
ENDIF
ENDDO
C
C ***** SELESAIKAN DENGAN TRIDIAGONAL SOLVER
C
CALL TRIDAG(A,B,C,D,U,IM-2)
DO I=1,IM-2
UB(I+1) = U(I)
V-14
ENDDO
IF( TIME.LT.TMAX )GOTO 100
C
C *** PRINT HASIL
C
WRITE(3,'(10X,999F10.3)') (TP(J),J=1,K)
DO I=1,IM
WRITE(3,'(1000F10.3)') XP(I),(UP(I,J),J=1,K)
ENDDO
STOP
END
C
C--------------------------< SUBROUTINE THOMAS >-------------------------------
C
SUBROUTINE tridag(a,b,c,r,u,n)
PARAMETER (NMAX=1000)
REAL a(n),b(n),c(n),r(n),u(n),bet,gam(NMAX)
IF( b(1).eq.0.) PAUSE 'tridag: rewrite equations'
c
c *** Forward Sweep
c
bet=b(1)
u(1)=r(1)/bet
do 11 j=2,n
gam(j)=c(j-1)/bet
bet=b(j)-a(j)*gam(j)
if(bet.eq.0.)pause 'tridag failed'
u(j)=(r(j)-a(j)*u(j-1))/bet
11 continue
c
c *** Backward Sweep
c
do 12 j=n-1,1,-1
u(j)=u(j)-gam(j+1)*u(j+1)
12 continue
RETURN
END

请先登录，再回答此问题。

回答（1 个）

Ben Barrowes
Ben Barrowes 2022-3-12
You should try my file exchange program, f2matlab. But you would have to convert it fortran90 style first, including refactoring the "goto" in there.
But I put it through my version of f2matlab after cleaning up the code a bit. I had to clean it up before it would compile, because when you copied and pasted this code, you lost all the indenting. And for example the ":" at the beginning of the line in your code is the continuation character. One thing I didn't get in this code was the "V-12", and "V-13", and "V-14". They almost look like page numbers. I commented them out.
Here is your converted code, which I admit I have not tested... nor can I without input files, etc. But there are support functions as well like writeFmt which I am including as an attachment.
function implisit(varargin)
global unit2fid;
if ~isempty(unit2fid);
unit2fid={};
end
imx = 1000;
a=zeros(1,imx);
b=zeros(1,imx);
c=zeros(1,imx);
cn=0.0;
cr=0.0;
d=zeros(1,imx);
dt=0.0;
dx=0.0;
i=0;
ii=0;
im=0;
ip=0;
iprint=0;
iter=0;
j=0;
k=0;
theta1=0.0;
theta2=0.0;
time_fv=0.0;
tmax=0.0;
tp=zeros(1,imx);
u=zeros(1,imx);
ub=zeros(1,imx);
ul=zeros(1,imx);
up=zeros(imx,100);
xl=0.0;
xp=zeros(1,imx);
xr=0.0;
% *** PROGRAM PENYELESAIAN PERSAMAAN HYPERBOLIK DENGAN METODA IMPLICIT
%------------------------------------------------------------------------------
%
% UL(I) = HARGA LAMA VARIABLE PADA TITIK I
% UB(I) = HARGA BARU VARIABLE PADA TITIK I
% V-12
% UP(I,K) = BEBERAPA HARGA VARIABLE PADA TITIK I DISIMPAN UNTUK PLOTTING
%
%------------------------------------------------------------------------------
%
if exist(strtrim('IM1.CON'))~=2;
fclose(fopen(strtrim('IM1.CON'),'w'));
end
thismlfid=fopen(strtrim('IM1.CON'),'r+');
unit2fid{end+1,1}=2;
unit2fid{end,2}=thismlfid;
unit2fid{end,3}=0;
unit2fid{end,4}=strtrim('IM1.CON');
unit2fid{end,5}=0;
unit2fid{end,6}=0;
% unit2fid maps fortran unit numbers to matlab fid numbers
if exist(strtrim('IM1.OUT'))~=2;
fclose(fopen(strtrim('IM1.OUT'),'w'));
end
thismlfid=fopen(strtrim('IM1.OUT'),'r+');
unit2fid{end+1,1}=3;
unit2fid{end,2}=thismlfid;
unit2fid{end,3}=0;
unit2fid{end,4}=strtrim('IM1.OUT');
unit2fid{end,5}=0;
unit2fid{end,6}=0;
% unit2fid maps fortran unit numbers to matlab fid numbers
%
% *** INPUT DATA DAN PARAMETER
%
im = 801;
xl = 0;
xr = 20;
dx =(xr-xl)./(im-1);
dt = 0.5;
tmax = 101;
cr = 0.1;
theta2 = 0.5;
theta1 = 1.0 - theta2;
iprint = 20;
cn = cr.*dt./dx;
%
% *** KOORDINAT NODE
%
for i = 1: im
xp(i) = xl +(i-1).*dx;
end
i = fix(im+1);
%
% *** HARGA AWAL (INITIAL VALUES)
%
for i = 1: im
if(xp(i) > 2.5 && xp(i) <= 3.0)
ub(i) =(xp(i)-2.5);
elseif(xp(i) > 3.0 && xp(i) <= 3.5)
ub(i) =(3.5-xp(i));
else;
ub(i) = 0.0;
end
end
i = fix(im+1);
%
% *** TETAPKAN HARGA UB(I) PADA SYARAT BATAS EXTERNAL
%
ub(1) = 0;
ub(im) = 0;
%
% *** LOOP ITERASI FTCS
%
iter = 0;
ip = fix(iprint);
time_fv = 0;
k = 0;
while (1);
% V-13
%
% ***** SIMPAN HARGA UB(I) PADA TIME STEP N
%
for i = 1: im
ul(i) = ub(i);
end
i = fix(im+1);
%
% **** SIMPAN HASIL UNTUK PLOTTING
%
if(ip == iprint)
ip = 0;
k = fix(k + 1);
for i = 1: im
up(i,k) = ub(i);
tp(k) = time_fv;
end
i = fix(im+1);
[writeErrFlag,wfso]=writeFmt(1,['%10i',repmat('%10.3f',1,3)],'iter','time_fv','cn');
eval(wfso);
[writeErrFlag,wfso]=writeFmt(2,['%10i',repmat('%10.3f',1,100)],'iter','time_fv',{'ub(i)','i',1,1,im});
eval(wfso);
end
iter = fix(iter + 1);
time_fv = time_fv + dt;
ip = fix(ip + 1);
%
% ***** PEMBENTUKAN MATRIK TRIDIAGONAL UNTUK MEMPERBAHARUI HARGA UB(I) PADA TIME STEP N+1
%
for i = 2: im - 1
ii = fix(i - 1);
if(i == 2)
a(ii) = 0;
b(ii) = 1.0;
c(ii) = theta1.*cn./2;
d(ii) = ul(i) + theta2.*cn./2.*(ul(i+1)-ul(i-1)) + theta1.*cn./2.*ul(i-1);
elseif(i == im-1)
a(ii) = -theta1.*cn./2;
b(ii) = 1.0;
c(ii) = 0;
d(ii) = ul(i) - theta2.*cn./2.*(ul(i+1)-ul(i-1)) - theta1.*cn./2.*ul(i+1);
else;
a(ii) = -theta1.*cn./2;
b(ii) = 1.0;
c(ii) = theta1.*cn./2;
d(ii) = ul(i) - theta2.*cn./2.*(ul(i+1)-ul(i-1));
end
end
i = fix(im - 1+1);
%
% ***** SELESAIKAN DENGAN TRIDIAGONAL SOLVER
%
[a,b,c,d,u,dumvar6]=tridag(a,b,c,d,u,im-2);
for i = 1: im - 2
ub(i+1) = u(i);
% V-14
end
i = fix(im - 2+1);
if(time_fv < tmax)
continue;
end
%
% *** PRINT HASIL
%
[writeErrFlag,wfso]=writeFmt(3,['%10x',repmat('%10.3f',1,999)],{'tp(j)','j',1,1,k});
eval(wfso);
for i = 1: im
[writeErrFlag,wfso]=writeFmt(3,[repmat('%10.3f',1,1000)],'xp(i)',{'up(i,j)','j',1,1,k});
eval(wfso);
end
i = fix(im+1);
tempBreak=1;
break;
end
closeAllOpenFiles(unit2fid);
end %program implisit

类别

在 Help Center 和 File Exchange 中查找有关 Fortran with MATLAB 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by