急救：关于多重打靶或并行打靶法的程序

贝贝 · 发表于 2006-7-16 09:56

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

我现在正在解一非线性刚性边值问题，查阅资料说要用多重打靶法（并行打靶法）来解，不知那位大侠有这样的程序，借小妹一用，感激不尽！
xiexiemx@163.com

[ 本帖最后由贝贝于 2006-7-16 10:46 编辑 ]

gghhjj · 发表于 2006-8-25 08:01

要什么语言的？matlab、fortran、c的我到是都收藏了一个，不过由于不是做这方面的，所以没测试

sczhang · 发表于 2006-8-25 16:40

呵呵，给我一份Ｃ的吧，我做混沌分叉方面的，我的邮箱：sczhang05@163.com,谢谢！！

gghhjj · 发表于 2006-8-26 01:14

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[1],yout[1];
int i,j,n;
double d,h=0.0;
double y[1],z[1],dydx[1],dzdx[1],yout[1],zout[1],yding[1],zding[1];
void derivs(double x,double y[],double z[],double dydx[],double dzdx[])
{
dydx[0]=z[0];//y'=z[0]
dzdx[0]=1.5*y[0]*y[0];//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[2],dzt[2],dzm[2],dzf[2],yt[2];
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[1],Fs[1],zjia[1],dlts=0.00001;
double Fpie[1],dltF;
n=1;
x=0;
y[0]=4;
z[0]=-1.0000;
//dydx[0]=z[0];//y'=z[0]
//dzdx[0]=1.5*y[0]*y[0];//z'=f(x,y,z)=1.5*y*y
do
{
for(j=0;j<5;j++)
{
if(j<1)
{
zding[0]=z[0];//zding[0]=z[0.0]
yding[0]=y[0];//yding[0]=y[0.0]
}
h=0.2+0.2*j;
rk4(n,h,x,y,z,dydx,dzdx,yout,zout);
////// cout<<"yout["<<h<<"]="<<yout[0]<<endl;//如果输出的yout[0]，也可以!!!!
}
F[0]=yout[0]-1.000000;//输出F(s);
zjia[0]=zding[0]+dlts;//zjia[0]=s(i)+dlts(i)
dydx[0]=zjia[0];//y'=z[0]
dzdx[0]=1.5*yding[0]*yding[0];//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[0]<<endl;
}
Fpie[0]=yout[0]-1.000000;
dltF=(Fpie[0]-F[0])/dlts;
z[0]=zding[0]-F[0]/dltF;//输出新的z[0] 也就是s(i+1)
Fs[0]=F[0];
//cout<<Fpie[0]<<endl;
//cout<<yout[0]<<endl;
// F[0]=Fpie[0];
y[0]=4;
cout<<"F[0]="<<F[0];
cout<<endl;
}while(F[0]>0.01);// &&F[0]<-0.00001);
cout<<"ok!"<<endl;
cout<<"z[0]="<<z[0]+Fs[0]/dltF<<endl;
cout<<"F[0]="<<F[0]<<endl;
cout<<"yout[0]="<<yout[0];
}

复制代码

转自：http://www.pudn.com/downloads29/sourcecode/math/detail93773.html

[ 本帖最后由 gghhjj 于 2007-1-8 06:34 编辑 ]

malong · 发表于 2006-8-26 08:35

能不能也把matlab编的也贴出来谢谢！我也想学习一下！

gghhjj · 发表于 2006-8-29 07:56

function dy=shootingfun(t,y);
% 定义打靶法的微分方程
% y''+ty'-4y=12t^2-3t 0<t<1
% y(0)=0,y(1)=2
dy=[y(2);4*y(1)-t.*y(2)+12*t.^2-3*t];

复制代码

% 关于线性打靶法的GUI文件,函数方程用shootingfun.m文件
close all
figure;
axes('position',[0.06,0.16,0.7,0.6]);
p='y''(0)=';P='\Delta(y(1))';Q=[];
Sl=uicontrol(gcf,'style','slider',...
'unit','normalized','position',[0.93,0.02,0.03,0.9],...
'BackgroundColor',[0.9 0.9 0.9],'ForegroundColor','r',...
'fontsize',14,'SliderStep',[0.01,0.01]);
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'',[p,q]);eval(SS);']);
Te=uicontrol(gcf,'style','text',...
'unit','normalized','position',[0.78,0.92,0.24,0.05],...
'BackgroundColor','w','ForegroundColor','r',...
'fontsize',10,'string',['y''(0)=',...
num2str(get(Sl,'value'))]);
Ed1=uicontrol(gcf,'style','edit',...
'unit','normalized','position',[0.81,0.86,0.1,0.05],...
'BackgroundColor','w','ForegroundColor','r',...
'string','2');
Ed0=uicontrol(gcf,'style','edit',...
'unit','normalized','position',[0.81,0.02,0.1,0.05],...
'BackgroundColor','w','ForegroundColor','r',...
'string','0');
Te0=uicontrol(gcf,'style','text',...
'unit','normalized','position',[0.11,0.8,0.4,0.2],...
'BackgroundColor',[0.8 0.8 0.9],'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([0,1],[0,2],'r*');hold on;plot([0,1],[0,2],'rs');
SS=['[t,y]=ode45(@shootingfun,[0,1],[0,y0]);',...
'set(h,''xdata'',t);set(h,''ydata'',y(:,1));',...
'set(DD,''string'',{P,num2str(y(end,1)-2)});'];
h=plot([0,1],[0,0]);
DD=uicontrol(gcf,'style','text',...
'unit','normalized','position',[0.55,0.8,0.18,0.14],...
'BackgroundColor','w','ForegroundColor','r',...
'string',{P,Q});

复制代码

转自：http://www.pudn.com/downloads29/sourcecode/math/detail93773.html

[ 本帖最后由 gghhjj 于 2007-1-8 06:35 编辑 ]

xjzuo · 发表于 2006-12-16 16:15

这个...转帖还是注明一下出处较好.

zqb138 · 发表于 2006-12-19 22:30

该matlab程序出现以下错误提示
??? Undefined function or variable 'Ed1'.

??? Error while evaluating uicontrol Callback.

shenyongjun · 发表于 2007-1-7 11:23

Matlab7上刚运行，没问题。谢谢楼主！

gghhjj · 发表于 2007-1-8 06:35

原帖由 xjzuo 于 2006-12-16 16:15 发表
这个...转帖还是注明一下出处较好.

谢谢您的意见，已经注明

gghhjj · 发表于 2007-1-8 06:37

原帖由 zqb138 于 2006-12-19 22:30 发表
该matlab程序出现以下错误提示
??? Undefined function or variable 'Ed1'.

??? Error while evaluating uicontrol Callback.

说明你的运行过程，应该是你使用问题

zjults · 发表于 2007-1-11 09:42

这些方法在gsl包里有，gnu science library.这些方法都可用于刚性常微分方程
如果新一点，可以用lie级数方法，李群积分方法。如果没有阻尼，可以用辛算法

shenfei · 发表于 2007-1-12 23:36

麻烦楼主给我一份fortran的吧，我做混沌分叉方面的，我的邮箱：shenf02@yahoo.com.cn,谢谢

gghhjj · 发表于 2007-1-14 01:58

!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 编辑 ]

gghhjj · 发表于 2007-1-14 01:59

原帖由 shenfei 于 2007-1-12 23:36 发表
麻烦楼主给我一份fortran的吧，我做混沌分叉方面的，我的邮箱：shenf02@yahoo.com.cn,谢谢

见14楼

账号		自动登录	找回密码
密码			我要加入

[分形与混沌] 急救：关于多重打靶或并行打靶法的程序

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

评分

评分

回复

关键是用embedded Runge-kutta 或 implict Runge-kutta方法

评分

评分