我有一個簡單的一階微分方程系統需要在Matlab中求解。它看起來像如下:Fortran |將函數作爲其他函數的參數傳遞
function passing
y0 = [0; 0];
t = 0:0.05:1;
y = myprocedure(@myfunc,0.05,t,y0);
myans = y'
end
function f = myfunc(t,y)
f = [-y(1) + y(2);
-0.8*t + y(2)];
end
function y=myprocedure(f,h,t,y0)
n = length(t);
y = zeros(length(y0),n);
y(:,1) = y0;
for i=1:n-1
k1 = f(t(i),y(:,i));
k2 = f(t(i)+h/2, y(:,i)+h/2*k1);
k3 = f(t(i)+h/2, y(:,i)+h/2*k2);
k4 = f(t(i)+h, y(:,i)+h*k3);
y(:,i+1) = y(:,i)+h/6*(k1+2*k2+2*k3+k4);
end
end
因此,運行passing.m文件後,我有結果:
>> passing
myans =
0 0
-0.0000 -0.0010
-0.0001 -0.0041
-0.0005 -0.0095
-0.0011 -0.0171
-0.0021 -0.0272
-0.0036 -0.0399
-0.0058 -0.0553
-0.0086 -0.0735
-0.-0.0946
-0.0169 -0.1190
-0.0225 -0.1466
-0.0293 -0.1777
-0.0374 -0.2124
-0.0469 -0.2510
-0.0579 -0.2936
-0.0705 -0.3404
-0.0849 -0.3917
-0.1012 -0.4477
-0.1196 -0.5086
-0.1402 -0.5746
現在,我想的passing.m代碼轉換爲Fortran 90的。我嘗試了很多東西,但我仍然迷失了方向 - 特別是在進行數組初始化並將函數傳遞給其他函數時。
任何幫助表示讚賞。謝謝!
==編輯:
我成功複製上面的Matlab代碼使用Fortran 90.我仍然不知道如何構建指針,但「接口」語句似乎工作。如果您查看代碼並提供一些輸入,我將不勝感激。謝謝。
module odemodule
implicit none
contains
function odef(a,b)
implicit none
real::a
real::b(2)
real::odef(2)
odef(1) = -b(1) + b(2)
odef(2) = -0.8*a + b(2)
end function odef
subroutine rk4(f,h,t,y0,y)
implicit none
real,dimension(2,21)::y
real,dimension(2)::k1,k2,k3,k4
real::h
real::t(21)
real::y0(2)
integer::i
interface
function f(a,b)
real::a
real::b(2), f(2)
end function f
end interface
y(1,1)=y0(1)
y(2,1)=y0(2)
do i=1,20
k1 = f(t(i),y(:,i))
k2 = f(t(i)+h/2, y(:,i)+h/2*k1)
k3 = f(t(i)+h/2, y(:,i)+h/2*k2)
k4 = f(t(i)+h, y(:,i)+h*k3)
y(:,i+1) = y(:,i)+h/6*(k1+2*k2+2*k3+k4)
end do
end subroutine rk4
end module odemodule
program passing
use odemodule
implicit none
real,parameter::dt=0.05
real::yinit(2)
real::y(2,21)
integer::i
real::t(21)
t(1)=0.0
do i=2,21
t(i)=t(i-1)+dt
end do
yinit(1)=0
yinit(2)=0
call rk4(odef,dt,t,yinit,y)
do i=1,21
write(*,100) y(1,i), y(2,i)
enddo
100 format(f7.4,3x,f7.4)
end program passing
@HighPerformanceMark - 沒問題。我花了很長時間才弄清楚了。 – mgilson