Программа для численного эксперимента 8.14

clear;

t0=0; tk=1;

a=1.5; b=1.5;

st=1.73;

nn=1;

for jj=3:7,

randn ('state', 0);

dt=1/(2^jj);

t=t0:dt:tk;

randn ('state', 0);

l=max(size(t))-1;

he=0;

for i=1:20*100,

x(1)=1; y(1)=1;

sw=0;

for p=1:l,

ksi=randn(1,nn+2);

sw=sw+ksi(1,1);

x(p+1)=x(1)*expm((a-0.5*(b^2))*p*dt+b*sqrt(dt)*sw);

ui0=sqrt(dt)*ksi(1,1);

ui1=-0.5*(dt^(3/2))*(ksi(1,1)+(1/sqrt(3))*ksi(1,2));

ui00=0.5*dt*((ksi(1,1))^2-1);

ui000=(1/6)*(dt^(3/2))*((ksi(1,1))^3-3*ksi(1,1));

ui0000=(1/24)*(dt^2)*((ksi(1,1))^4-6*((ksi(1,1))^2)+3);

bbb=(2/3)*((ksi(1,1))^2)+(1/sqrt(3))*ksi(1,1)*ksi(1,2)-(1/(3*sqrt(5)))*ksi(1,3)*ksi(1,1)-1;

for ij=1:nn,

bbb=bbb+((ksi(1,ij))^2)/((2*ij-1)*(2*ij+3))-(ksi(1,ij)*ksi(1,ij+2))/((2*ij+3)*sqrt((2*ij+1)*(2*ij+5)));

end;

ui10=-(1/4)*(dt^2)*bbb;

ui01=0.5*(dt^2)+ui1*ui0-ui10;

z=y(p)+sqrt(dt)*b*y(p);

u=y(p)-sqrt(dt)*b*y(p);

v=y(p)-3*sqrt(dt)*b*y(p);

w=y(p)+3*sqrt(dt)*b*y(p);

r=y(p)+dt*a*y(p);

f=y(p)+2*dt*b*y(p);

g=y(p)+0.5*dt*a*y(p);

q=y(p)-0.5*dt*a*y(p);

uu=b*y(p)*ui0+(1/sqrt(dt))*((1/48)*b*v-(1/48)*b*w+(27/48)*b*z-(27/48)*b*u)*ui00;

vv=(1/dt)*((0.5*a*f-0.5*a*y(p))*(0.5*dt*ui0+ui1)-(0.5*(b*z+b*u-2*b*y(p))+b*g-b*q)*ui1);

z11=b*(v+0.5*sqrt(dt)*b*v)-b*(v-0.5*sqrt(dt)*b*v);

z12=b*(w+0.5*sqrt(dt)*b*w)-b*(w-0.5*sqrt(dt)*b*w);

z21=b*(z+0.5*sqrt(dt)*b*z)-b*(z-0.5*sqrt(dt)*b*z);

z22=b*(u+0.5*sqrt(dt)*b*u)-b*(u-0.5*sqrt(dt)*b*u);

zzz=(1/dt)*((1/48)*z11-(1/48)*z12+(27/48)*z21-(27/48)*z22)*ui000;

u11=0.5*(b*v+b*(v+2*sqrt(dt)*b*v)-2*b*(v+sqrt(dt)*b*v))+b*(v+0.5*dt*a*v)-b*(v-0.5*dt*a*v);

u12=0.5*(b*w+b*(w+2*sqrt(dt)*b*w)-2*b*(w+sqrt(dt)*b*w))+b*(w+0.5*dt*a*w)-b*(w-0.5*dt*a*w);

u21=0.5*(b*z+b*(z+2*sqrt(dt)*b*z)-2*b*(z+sqrt(dt)*b*z))+b*(z+0.5*dt*a*z)-b*(z-0.5*dt*a*z);

u22=0.5*(b*u+b*(u+2*sqrt(dt)*b*u)-2*b*(u+sqrt(dt)*b*u))+b*(u+0.5*dt*a*u)-b*(u-0.5*dt*a*u);

uuu=((1/48)*u11-(1/48)*u12+(27/48)*u21-(27/48)*u22)*(ui10-ui01);

h11=0.5*a*(v+2*dt*b*v)-0.5*a*v;

h12=0.5*a*(w+2*dt*b*w)-0.5*a*w;

h21=0.5*a*(z+2*dt*b*z)-0.5*a*z;

h22=0.5*a*(u+2*dt*b*u)-0.5*a*u;

www=((1/48)*h11-(1/48)*h12+(27/48)*h21-(27/48)*h22)*(ui01+0.5*dt*ui00);

uk11=0.5*b*(z+2*sqrt(dt)*b*z)-0.5*b*z;

uk12=0.5*b*(u+2*sqrt(dt)*b*u)-0.5*b*u;

uk21=0.5*b*(y(p)+2*sqrt(dt)*b*y(p))-0.5*b*y(p);

uk22=0.5*b*(g+2*sqrt(dt)*b*g)-0.5*b*g;

uk23=0.5*b*(q+2*sqrt(dt)*b*q)-0.5*b*q;

sss=-(0.5*(uk11+uk12-2*uk21)+uk22-uk23)*ui10;

r11=v+0.5*sqrt(dt)*b*v;

r12=v-0.5*sqrt(dt)*b*v;

f11=w+0.5*sqrt(dt)*b*w;

f12=w-0.5*sqrt(dt)*b*w;

g11=z+0.5*sqrt(dt)*b*z;

g12=z-0.5*sqrt(dt)*b*z;

q11=u+0.5*sqrt(dt)*b*u;

q12=u-0.5*sqrt(dt)*b*u;

c11=0.5*b*(r11+2*sqrt(dt)*b*r11)-0.5*b*r11-0.5*b*(r12+2*sqrt(dt)*b*r12)+0.5*b*r12;

c12=0.5*b*(f11+2*sqrt(dt)*b*f11)-0.5*b*f11-0.5*b*(f12+2*sqrt(dt)*b*f12)+0.5*b*f12;

c21=0.5*b*(g11+2*sqrt(dt)*b*g11)-0.5*b*g11-0.5*b*(g12+2*sqrt(dt)*b*g12)+0.5*b*g12;

c22=0.5*b*(q11+2*sqrt(dt)*b*q11)-0.5*b*q11-0.5*b*(q12+2*sqrt(dt)*b*q12)+0.5*b*q12;

ccc=((1/48)*c11-(1/48)*c12+(27/48)*c21-(27/48)*c22)*ui0000;

y(p+1)=(1/(1-0.5*a*dt))*(y(p)+0.5*a*y(p)*dt+uu+vv+zzz+(uuu+www+sss+ccc)*(1/(dt^(3/2))));

end;

e(i)=abs(x(l+1)-y(l+1));

he=he+e(i);

end;

he=he/(20*100);

for j=1:20,

ee=0;

for k=(j-1)*100+1:j*100,

ee=ee+e(k);

end;

hh(j)=ee/100;

end;

s=0;

for j=1:20,

s=s+(hh(j)-he)^2;

end;

s=s/19;

hee(jj-2)=he;

he1(jj-2)=he-st*sqrt(s/20);

he2(jj-2)=he+st*sqrt(s/20);

ttt(jj-2)=dt;

end;

e1=[he1(1) he2(1)];

e2=[he1(2) he2(2)];

e3=[he1(3) he2(3)];

e4=[he1(4) he2(4)];

e5=[he1(5) he2(5)];

tt1=[1/(2^3) 1/(2^3)];

tt2=[1/(2^4) 1/(2^4)];

tt3=[1/(2^5) 1/(2^5)];

tt4=[1/(2^6) 1/(2^6)];

tt5=[1/(2^7) 1/(2^7)];

plot(tt1,e1,tt2,e2,tt3,e3,tt4,e4,tt5,e5,ttt,hee);

pause;

Рис. 17.4 Результат численного эксперимента 8.14