急救：关于多重打靶或并行打靶法的程序
我现在正在解一非线性刚性边值问题，查阅资料说要用多重打靶法（并行打靶法）来解，不知那位大侠有这样的程序，借小妹一用，感激不尽！xiexiemx@163.com
[ 本帖最后由 贝贝 于 2006-7-16 10:46 编辑 ] 要什么语言的？matlab、fortran、c的我到是都收藏了一个，不过由于不是做这方面的，所以没测试 呵呵，给我一份Ｃ的吧，我做混沌分叉方面的，我的邮箱：sczhang05@163.com,谢谢！！ C语言写的打靶法程序，用于数值计算中的边值问题，程序中采用rugga-kutta法对常微分方程处理（程序未经测试，如有问题请见谅）#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>
# include<malloc.h>
# include<stdio.h>
# include<float.h>
double x;
double hh,h6,xh;//,dydx,yout;
int i,j,n;
double d,h=0.0;
double y,z,dydx,dzdx,yout,zout,yding,zding;
void derivs(double x,double y[],double z[],double dydx[],double dzdx[])
{
dydx=z;//y'=z
dzdx=1.5*y*y;//z'=f(x,y,z)=1.5*y*y
}
//子过程RK4
void rk4(int n,double h,double x,double y[],double z[],double dydx[],double dzdx[],double yout[],double zout[])
{
double zt,dzt,dzm,dzf,yt;
hh=h*0.5;
h6=h/6;
xh=x+hh;
derivs(x,y,z,dydx,dzdx);
for(i=0;i<n;i++)
{
zt=z+hh*dzdx;//z(n)+h/2*L1
yt=y+hh*z;
}
derivs(xh,yt,zt,dydx,dzt);//x+h/2,y(n)+h/2*zn,zn+h/2*L1,L2
for(i=0;i<n;i++)
{
yt=y+hh*z+h*h/4*dzdx;
zt=z+hh*dzt;//z(n)+h/2*L2
}
derivs(xh,yt,zt,dydx,dzm);//dzm=L3
for(i=0;i<n;i++)
{
yt=y+h*z+h*h/2*dzt;//dzt=L2
zt=z+h*dzm;//z(n)+h/2*L3
//dym=dyt+dym;///K2+k3
}
derivs(x+h,yt,zt,dydx,dzf);/////x+h*******dzf=L4
for(i=0;i<n;i++)
{
yout=y+h*z+h*h/6*(dzdx+dzt+dzm);
y=yout;
zout=z+h6*(dzdx+2.0*dzt+2.0*dzm+dzf); //h/6=h6
// cout<<zout<<endl;//yout从这里输出,切记！！！// cout<<h6;
z=zout;
}
for(i=0;i<n;i++)
{
dzm=0.0;
dzt=0.0;
yt=0.0;
zt=0.0;
dzf=0.0;
}
}
void main()
{
double F,Fs,zjia,dlts=0.00001;
double Fpie,dltF;
n=1;
x=0;
y=4;
z=-1.0000;
//dydx=z;//y'=z
//dzdx=1.5*y*y;//z'=f(x,y,z)=1.5*y*y
do
{
for(j=0;j<5;j++)
{
if(j<1)
{
zding=z;//zding=z
yding=y;//yding=y
}
h=0.2+0.2*j;
rk4(n,h,x,y,z,dydx,dzdx,yout,zout);
////// cout<<"yout["<<h<<"]="<<yout<<endl;//如果输出的yout，也可以!!!!
}
F=yout-1.000000;//输出F(s);
zjia=zding+dlts;//zjia=s(i)+dlts(i)
dydx=zjia;//y'=z
dzdx=1.5*yding*yding;//z'=f(x,y,z)=1.5*y*y
for(j=0;j<5;j++)
{
h=0.2+0.2*j;
rk4(n,h,x,yding,zjia,dydx,dzdx,yout,zout);
// cout<<"yout["<<h<<"]="<<yout<<endl;
}
Fpie=yout-1.000000;
dltF=(Fpie-F)/dlts;
z=zding-F/dltF;//输出新的z 也就是s(i+1)
Fs=F;
//cout<<Fpie<<endl;
//cout<<yout<<endl;
// F=Fpie;
y=4;
cout<<"F="<<F;
cout<<endl;
}while(F>0.01);// &&F<-0.00001);
cout<<"ok!"<<endl;
cout<<"z="<<z+Fs/dltF<<endl;
cout<<"F="<<F<<endl;
cout<<"yout="<<yout;
}
转自：http://www.pudn.com/downloads29/sourcecode/math/detail93773.html
[ 本帖最后由 gghhjj 于 2007-1-8 06:34 编辑 ] 能不能也把matlab编的 也贴出来谢谢！ 我也想学习一下！ function dy=shootingfun(t,y);
% 定义打靶法的微分方程
% y''+ty'-4y=12t^2-3t 0<t<1
% y(0)=0,y(1)=2
dy=;
% 关于线性打靶法的GUI文件,函数方程用shootingfun.m文件
close all
figure;
axes('position',);
p='y''(0)=';P='\Delta(y(1))';Q=[];
Sl=uicontrol(gcf,'style','slider',...
'unit','normalized','position',,...
'BackgroundColor',,'ForegroundColor','r',...
'fontsize',14,'SliderStep',);
set(Sl,'callback',['a=str2num(get(Ed1,''string''));',...
'b=str2num(get(Ed0,''string''));',...
'y0=b+(a-b)*get(Sl,''value'');q=num2str(y0);',...
'set(Te,''string'',);eval(SS);']);
Te=uicontrol(gcf,'style','text',...
'unit','normalized','position',,...
'BackgroundColor','w','ForegroundColor','r',...
'fontsize',10,'string',['y''(0)=',...
num2str(get(Sl,'value'))]);
Ed1=uicontrol(gcf,'style','edit',...
'unit','normalized','position',,...
'BackgroundColor','w','ForegroundColor','r',...
'string','2');
Ed0=uicontrol(gcf,'style','edit',...
'unit','normalized','position',,...
'BackgroundColor','w','ForegroundColor','r',...
'string','0');
Te0=uicontrol(gcf,'style','text',...
'unit','normalized','position',,...
'BackgroundColor',,'ForegroundColor','b',...
'string',{'differential equations:',...
'y''''+ty''-4y=12t^2-3t 0<t<1','y(0)=0,y(1)=2'},...
'fontsize',14,'HorizontalAlignment','left');
plot(,,'r*');hold on;plot(,,'rs');
SS=['=ode45(@shootingfun,,);',...
'set(h,''xdata'',t);set(h,''ydata'',y(:,1));',...
'set(DD,''string'',{P,num2str(y(end,1)-2)});'];
h=plot(,);
DD=uicontrol(gcf,'style','text',...
'unit','normalized','position',,...
'BackgroundColor','w','ForegroundColor','r',...
'string',{P,Q});
转自：http://www.pudn.com/downloads29/sourcecode/math/detail93773.html
[ 本帖最后由 gghhjj 于 2007-1-8 06:35 编辑 ]
回复
这个...转帖还是注明一下出处较好. 该matlab程序出现以下错误提示??? Undefined function or variable 'Ed1'.
??? Error while evaluating uicontrol Callback. Matlab7上刚运行，没问题。谢谢楼主！ 原帖由 xjzuo 于 2006-12-16 16:15 发表
这个...转帖还是注明一下出处较好.
谢谢您的意见，已经注明 原帖由 zqb138 于 2006-12-19 22:30 发表
该matlab程序出现以下错误提示
??? Undefined function or variable 'Ed1'.
??? Error while evaluating uicontrol Callback.
说明你的运行过程，应该是你使用问题
关键是用embedded Runge-kutta 或 implict Runge-kutta方法
这些方法在gsl包里有，gnu science library.这些方法都可用于刚性常微分方程如果新一点，可以用lie级数方法，李群积分方法。如果没有阻尼，可以用辛算法 麻烦楼主给我一份fortran的吧，我做混沌分叉方面的，我的邮箱：shenf02@yahoo.com.cn,谢谢 !C *********************************************************
!C * fixed-fixed *
!C *********************************************************
PARAMETER (NVAR=8,N2=3,DELTA=0.1E-1,EPS=0.1E-3)
DIMENSION V(3),DELV(3),F(3),DV(3),YS(8)
COMMON HL,Wmax,T
OPEN(1,FILE='xb1.DAT')
OPEN(2,FILE='xb2.DAT')
OPEN(3,FILE='xb3.DAT')
OPEN(4,FILE='xb4.DAT')
OPEN(5,FILE='xb5.DAT')
!WRITE(*,*) ' V(1),V(2),V(3)='
!READ(*,*)V(1),V(2),V(3)
V(1)=0.1
V(2)=0.2
V(3)=1.0
HL=40.0
T=50.0
X1=0.5
H1=0.1
HMIN=1.0E-30
X2=1.0
DO 10 J=10,150,1
W=J/100.0
Wmax=W/HL
DELV(1)=DELTA*V(1)
DELV(2)=DELTA*V(2)
DELV(3)=DELTA*V(3)
9 CALL SHOOT(NVAR,YS,V,DELV,N2,X1,X2,EPS,H1,HMIN,F,DV)
WRITE(*,*) W,V(1),V(2),v(3)
IF(ABS(DV(1)).GT.ABS(EPS*V(1)).OR.ABS(DV(2)).GT.ABS(EPS*V(2)).OR.ABS(DV(3)).GT.ABS(EPS*V(3))) GOTO 9
WRITE(1,*) W,V(1),V(2),v(3)
WRITE(2,*) V(3),W
WRITE(3,*) v(3),v(1)
WRITE(5,*) v(3),V(2)
10 CONTINUE
END
SUBROUTINE LOAD(X1,V,Y)
DIMENSION V(3),Y(8)
COMMON HL,Wmax,T
Y(1)=0.0
Y(2)=0.0
Y(3)=Wmax
Y(4)=0.0
Y(5)=V(1)
Y(6)=V(2)
Y(7)=0.0
Y(8)=V(3)
RETURN
END
SUBROUTINE SCORE(X2,Y,F)
DIMENSION Y(8),F(3)
COMMON HL,Wmax,T
F(1)=Y(2)
F(2)=Y(3)
F(3)=Y(4)
RETURN
END
!C
SUBROUTINE DERIVS(X,Y,DYDX)
DIMENSION Y(8),DYDX(8)
COMMON HL,Wmax,T
R=(T-Y(6)*COS(Y(4))-Y(7)*SIN(Y(4)))/(12.0*HL**2)+1.0
DYDX(1)=R
DYDX(2)=R*COS(Y(4))-1.0
DYDX(3)=R*SIN(Y(4))
DYDX(4)=Y(5)
DYDX(5)=R*(-Y(6)*sin(Y(4))+Y(7)*cos(Y(4)))
DYDX(6)=0.0
DYDX(7)=R*Y(8)
DYDX(8)=0.0
RETURN
END
SUBROUTINE SHOOT(NVAR,YS,V,DELV,N2,X1,X2,EPS,H1,HMIN,F,DV)
EXTERNAL DERIVS,RKQC
PARAMETER (NP=20)
DIMENSION V(N2),DELV(N2),F(N2),DV(N2),Y(NP),DFDV(NP,NP),INDX(NP),YS(NVAR)
CALL LOAD(X1,V,Y)
CALL ODEINT(Y,NVAR,X1,X2,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
CALL SCORE(X2,Y,F)
DO 12 IV=1,N2
SAV=V(IV)
V(IV)=V(IV)+DELV(IV)
CALL LOAD(X1,V,Y)
CALL ODEINT(Y,NVAR,X1,X2,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
DO 222 J=1,NVAR
222 YS(J)=Y(J)
CALL SCORE(X2,Y,DV)
DO 11 I=1,N2
DFDV(I,IV)=(DV(I)-F(I))/DELV(IV)
11 CONTINUE
V(IV)=SAV
12 CONTINUE
DO 13 IV=1,N2
DV(IV)=-F(IV)
13 CONTINUE
CALL LUDCMP(DFDV,N2,NP,INDX,DET)
CALL LUBKSB(DFDV,N2,NP,INDX,DV)
DO 14 IV=1,N2
V(IV)=V(IV)+DV(IV)
14 CONTINUE
RETURN
END
SUBROUTINE ODEINT(YSTART,NVAR,X1,X2,EPS,H1,HMIN,NOK,NBAD,DERIVS,RKQC)
PARAMETER (MAXSTP=10000,NMAX=20,TWO=2.0,ZERO=0.0,TINY=1.E-30)
COMMON /PATH/KMAX,KOUNT,DXSAV,XP(350),YP(15,350)
DIMENSION YSTART(NVAR),YSCAL(NMAX),Y(NMAX),DYDX(NMAX)
X=X1
H=SIGN(H1,X2-X1)
NOK=0
NBAD=0
KOUNT=0
DO 11 I=1,NVAR
Y(I)=YSTART(I)
11 CONTINUE
XSAV=X-DXSAV*TWO
DO 16 NSTP=1,MAXSTP
CALL DERIVS(X,Y,DYDX)
DO 12 I=1,NVAR
YSCAL(I)=ABS(Y(I))+ABS(H*DYDX(I))+TINY
12 CONTINUE
IF (KMAX.GT.0) THEN
IF (ABS(X-XSAV).GT.ABS(DXSAV)) THEN
IF (KOUNT.LT.KMAX-1) THEN
KOUNT=KOUNT+1
XP(KOUNT)=X
DO 13 I=1,NVAR
YP(I,KOUNT)=Y(I)
13 CONTINUE
XSAV=X
ENDIF
ENDIF
ENDIF
IF ((X+H-X2)*(X+H-X1).GT.ZERO) H=X2-X
CALL RKQC(Y,DYDX,NVAR,X,H,EPS,YSCAL,HDID,HNEXT,DERIVS)
IF (HDID.EQ.H) THEN
NOK=NOK+1
ELSE
NBAD=NBAD+1
ENDIF
IF ((X-X2)*(X2-X1).GE.ZERO) THEN
DO 14 I=1,NVAR
YSTART(I)=Y(I)
14 CONTINUE
IF (KMAX.NE.0) THEN
KOUNT=KOUNT+1
XP(KOUNT)=X
DO 15 I=1,NVAR
YP(I,KOUNT)=Y(I)
IF(I.GT.1) GO TO 15
15 CONTINUE
ENDIF
RETURN
ENDIF
IF (ABS(HNEXT).LT.HMIN) PAUSE 'Stepsize smaller then minimun.'
H=HNEXT
16 CONTINUE
PAUSE'Too many steps.'
RETURN
END
SUBROUTINE RKQC(Y,DYDX,N,X,HTRY,EPS,YSCAL,HDID,HNEXT,DERIVS)
PARAMETER (NMAX=20,ONE=1.,SAFETY=0.9,ERRCON=6.E-4, FCOR=6.6667E-2)
EXTERNAL DERIVS
DIMENSION Y(N),DYDX(N),YSCAL(N),YTEMP(NMAX),YSAV(NMAX),DYSAV(NMAX)
PGROW=-0.20
PSHRNK=-0.25
XSAV=X
DO 11 I=1,N
YSAV(I)=Y(I)
DYSAV(I)=DYDX(I)
11 CONTINUE
H=HTRY
1 HH=0.5*H
CALL RK4(YSAV,DYSAV,N,XSAV,HH,YTEMP,DERIVS)
X=XSAV+HH
CALL DERIVS(X,YTEMP,DYDX)
CALL RK4(YTEMP,DYDX,N,X,HH,Y,DERIVS)
X=XSAV+H
IF (X.EQ.XSAV) PAUSE 'Stepsize not significant in RKQC.'
CALL RK4(YSAV,DYSAV,N,XSAV,H,YTEMP,DERIVS)
ERRMAX=0.
DO 12 I=1,N
YTEMP(I)=Y(I)-YTEMP(I)
ERRMAX=MAX(ERRMAX,ABS(YTEMP(I)/YSCAL(I)))
12 CONTINUE
ERRMAX=ERRMAX/EPS
IF (ERRMAX.GT.ONE) THEN
H=SAFETY*H*(ERRMAX**PSHRNK)
GOTO 1
ELSE
HDID=H
IF (ERRMAX.GT.ERRCON) THEN
HNEXT=SAFETY*H*(ERRMAX**PGROW)
ELSE
HNEXT=4.*H
ENDIF
ENDIF
DO 13 I=1,N
Y(I)=Y(I)+YTEMP(I)*FCOR
13 CONTINUE
!write(4,*) X,y(1)
!write(5,*) 1-X,y(1)
RETURN
END
SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,DERIVS)
PARAMETER (NMAX=20)
DIMENSION Y(N),DYDX(N),YOUT(N),YT(NMAX),DYT(NMAX),DYM(NMAX)
HH=H*0.5
H6=H/6.
XH=X+HH
DO 11 I=1,N
YT(I)=Y(I)+HH*DYDX(I)
11 CONTINUE
CALL DERIVS(XH,YT,DYT)
DO 12 I=1,N
YT(I)=Y(I)+HH*DYT(I)
12 CONTINUE
CALL DERIVS(XH,YT,DYM)
DO 13 I=1,N
YT(I)=Y(I)+H*DYM(I)
DYM(I)=DYT(I)+DYM(I)
13 CONTINUE
CALL DERIVS(X+H,YT,DYT)
DO 14 I=1,N
YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I))
14 CONTINUE
RETURN
END
SUBROUTINE LUBKSB(A,N,NP,INDX,B)
DIMENSION A(NP,NP),INDX(N),B(N)
II=0
DO 12 I=1,N
LL=INDX(I)
SUM=B(LL)
B(LL)=B(I)
IF (II.NE.0) THEN
DO 11 J=II,I-1
SUM=SUM-A(I,J)*B(J)
11 CONTINUE
ELSE IF (SUM.NE.0.) THEN
II=I
ENDIF
B(I)=SUM
12 CONTINUE
DO 14 I=N,1,-1
SUM=B(I)
IF (I.LT.N) THEN
DO 13 J=I+1,N
SUM=SUM-A(I,J)*B(J)
13 CONTINUE
ENDIF
B(I)=SUM/A(I,I)
14 CONTINUE
RETURN
END
SUBROUTINE LUDCMP(A,N,NP,INDX,D)
PARAMETER (NMAX=100,TINY=1.0E-20)
DIMENSION A(NP,NP),INDX(N),VV(NMAX)
D=1.
DO 12 I=1,N
AAMAX=0.
DO 11 J=1,N
IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
11 CONTINUE
IF (AAMAX.EQ.0.) PAUSE'Singular matrix.'
VV(I)=1./AAMAX
12 CONTINUE
DO 19 J=1,N
IF (J.GT.1) THEN
DO 14 I=1,J-1
SUM=A(I,J)
IF (I.GT.1) THEN
DO 13 K=1,I-1
SUM=SUM-A(I,K)*A(K,J)
13 CONTINUE
A(I,J)=SUM
ENDIF
14 CONTINUE
ENDIF
AAMAX=0.
DO 16 I=J,N
SUM=A(I,J)
IF (J.GT.1) THEN
DO 15 K=1,J-1
SUM=SUM-A(I,K)*A(K,J)
15 CONTINUE
A(I,J)=SUM
ENDIF
DUM=VV(I)*ABS(SUM)
IF (DUM.GE.AAMAX) THEN
IMAX=I
AAMAX=DUM
ENDIF
16 CONTINUE
IF (J.NE.IMAX) THEN
DO 17 K=1,N
DUM=A(IMAX,K)
A(IMAX,K)=A(J,K)
A(J,K)=DUM
17 CONTINUE
D=-D
VV(IMAX)=VV(J)
ENDIF
INDX(J)=IMAX
IF (J.NE.N) THEN
IF (A(J,J).EQ.0.) A(J,J)=TINY
DUM=1./A(J,J)
DO 18 I=J+1,N
A(I,J)=A(I,J)*DUM
18 CONTINUE
ENDIF
19 CONTINUE
IF (A(N,N).EQ.0.) A(N,N)=TINY
RETURN
END
上述程序均来自pudn程序员联合开发网
[ 本帖最后由 gghhjj 于 2007-1-14 02:01 编辑 ] 原帖由 shenfei 于 2007-1-12 23:36 发表
麻烦楼主给我一份fortran的吧，我做混沌分叉方面的，我的邮箱：shenf02@yahoo.com.cn,谢谢
见14楼
页:
[1]
2