(*second order IDE*)
(*data*)
jan[x_,s_]:={{If[s<=x,(1-x)s,x(1-s)]}};
g[x_,u_]:={Exp[-u[[1]]]}; lamb={1};
(*module*)
dane=solveIDE[g,1,0,1,9,bound->{{0,1}}]
(*Qh(x) module*)
intPol[n_,p_]:=Module[{r,pary,vek};
Clear[x];
pary = Table[Table[{dane[[i,1]],dane[[i,r]]}, {i,1,n+2}], {r,2,p+1}];
vek[r_]:=Part[pary,r];
Table[Expand[InterpolatingPolynomial[vek[k],x]],
{k,1,p}]
];
(*command*)
intPol[9,1]
(*residual error*)
resError[g_,p_,n_]:=Module[{p1,pxx,ps,ca1,prawa};
p1[x_]:=intPol[n,p];
pxx=D[D[intPol[n,p],x],x];
ps=intPol[n,p]/.x- >s;
ca1=Expand[Integrate[jan[x,s].ps*lamb,{s,0,1}]];
prawa=g[x,p1[x]];
blad=-pxx+ca1-prawa;
Table[N[Prepend[blad,x]],{x,0.1,0.9,0.1} ]//TableForm
];
(*a system of IDE*)
jan[x_,s_]:={{s*sin[x],2*Exp[-s]}; {s,cos[x]*Exp[-s]}};
lamb = {1,1};
g[x_,u_]:={-1 u[[1]]-Exp[-x],-1 u[[2]]-4};
N[dane=solveIDE[g,2,0,1,9,bound- >{{0,1}; {0,1}}]]
(*command*)
intPol[9,2]
(*module solveIDE*)
Options[solveIDE] = { bound -> homogeneous,iters -> twon};
solveIDE[g_, p_, a_, b_, n_, opts_] := Module[
{ h, x, boundary,w0,v0, iterNumber, beta,
gamma,solution, bcorrections, useboundary,solve3, solve5, oneStep},
h = (b-a)/(n+1);
x = Table [ N[a + i h], {i, 1, n} ];
homogeneous = Table [0, {p}, {2} ];
twon = 2n;
boundary=bound/.{opts} /.Options[solveIDE];
v0=Table [(boundary[[q,2]]-boundary[[q,1]])*
i*h/(b-a)+boundary[[q,1]],{q,1,p},{i,1,n}];
iterNumber = iters/.{opts} /.Options[solveIDE];
beta={0};
Do[AppendTo[beta, 1/(14-Last[beta]) ], {3} ];
Do[AppendTo[beta, N[7 - 4 Sqrt[3]] ], {n-3} ];
gamma=Table [i/(i+1), {i, 1, n} ];
bcorrections=Table [
{{12*h^2* boundary[[q,1]],
- h^2*boundary[[q,1]]},{-h^2*boundary[[q,2]],
12*h^2*boundary[[q,2]]} } /h^2, {q,1,p} ];
useboundary[f_, bcorr ]:=Join[Take[f,2] +
bcorr[[1]],Take[f, {3,-3} ],
Take[f,-2] + bcorr[[2]]];
solve3[f_, d11_, alpha_]:= Module[{ f1, sol },
f1[1]=f[[1]]=d11;
f1[i_]:=f1[i]=(f[[i]]+f1[i-1])*alpha[[i]];
sol[n]=If[d11==12, f[[n]]/12, f1[n]];
sol[i_]:=sol[i]=f1[i]+alpha[[i]]*sol[i+1];
Table [sol[i],{i,1,n} ]];
solve5[f_]:=With[{w= solve3[f, 12, beta]},
solve3[w, 2, gamma]];
cab[i_,q_,r_,vv_]:=Module[{ks0,ks1,ks2,ma },
ks0=-Sqrt [3/5]; ks1=0; ks2=Sqrt[3/5];
N[Sum[5/9 jan[a+i h,a+j h+(1+ks0) h/2][[q,r]]*
((1-ks0)*vv[[r,j+1]]+(1+ks0)vv[[r,j+2]])+
8/9 jan[a+i h,a+j h+(1-ks1) h/2][[q,r]] *
((1-ks1)*vv[[r,j+1]]+(1+ks1)vv[[r,j+2]])+
5/9 jan[a+i h,a+j h+(1+ks2)h/2][[q,r]]*
((1-ks2)*vv[[r,j+1]]+(1+ks2)vv[[r,j+2]]),
{j,0,n} ]*h/4,4]
];
calka[ ii_,qq_,ww_]:=Sum[cab[ii,qq,rr,ww],{rr,1,p} ];
gaussint[i_,q_,r_,vv_]:=Module[{ },
m3=Mod[n+1,3];
last4=If[m3==0, n-1,n-m3-1];
integrandValues[k_, y_]:=Module[
{intpoints, intpol, integrand, c=Sqrt [3/5]},
intpoints = Table [{a+(j-1)*h,y[[j]]}, {j,k,k+3} ];
intpoly = InterpolatingPolynomial[intpoints,s] ;
integrand=jan[a+i*h, s][[q,r]]*intpoly;
With[ {xj=a+(j-1)*h},
Table [N[{ (5/9)*integrand/.s->xj+(1-c)*h/2,
(8/9)*integrand/.s->xj+ h/2,
(5/9)*integrand/.s->xj+(1+c)* h/2 } ],
{j, k,k+2} ]* h/2 ]];
integral3[k_,y_]:=Apply[Plus,
Flatten[integrandValues[k,y]]];
int1[y_]:=Apply[Plus, Map[integral3[#,y]&,
Range[1, last4,3]]];
restint[y_]:=If[m3 == 0, 0,
Apply[Plus, Flatten[
Take[integrandValues[n-1,y],-m3]]]
];
y=vv[[r]];
int1[y]+restint[y]];
kcalka[i_,q_,ww_]:=Sum[gaussint[i,q,r,ww],{r,1,p} ];
oneStep[v_]:=Module[ {ff, bff,d,q,r,s,y,zz},
ff=12*h^2*N[Transpose[
MapThread[g,{x,Transpose[v]} ]]];
zz=Transpose[boundary];
cc=Table [Append[Prepend[v[[q]],
zz[[1,q]]],zz[[2,q]]],{q,1,p} ];
fff=Table [lamb[[q]]*kcalka[i,q,cc],{i,1,n},{q,1,p} ];
ff=ff-12 h^2 Transpose[fff];
bff=MapThread[useboundary,{ff,bcorrections} ];
Map[solve5, bff] ];
solution = Nest[oneStep, v0, iterNumber];
PrependTo[solution, x];
With[ { initiala} =
Prepend[Table [boundary[[q,1]],
{q,1,p} ],a], endb=
Prepend[Table [boundary[[q,2]], {q,1,p} ],b],
Append[Prepend[Transpose[solution], initiala],endb]]
]
评论2