Mathematica代码_mathematica_有限差分法_

(*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]] ]