I am having error in this programme. i gave the call function also.

1 visualización (últimos 30 días)
GEETHA NANIAPPAN
GEETHA NANIAPPAN el 15 de Sept. de 2021
Comentada: Walter Roberson el 15 de Sept. de 2021
Nx=41;XL=0;XR=2*pi;
dx=(XR-XL)/(Nx-1);
for i=1:Nx
XX(i)=XL+(i-1)*dx;
end
cfl=0.125;dt=cfl*dx^3;tlast=1;
Nstep=tlast/dt;
cc=1;
IRK=4;
[u0ex]=feval(@JSC1,XX,0,cc);
uold=u0ex;
subplot(2,1,1),plot(XX,u0ex,'k-');
hold on
for k=1:Nstep
if IRK==4
u0x=reconuxxxp(uold,Nx,dx);
k0=dt*resJSC1(uold,u0x,cc);
u1=uold+k0/2;
u1X=reconuxxxp(u1,Nx,dx);
k1=dt*resJSC1(uold,u1x,cc);
u2=uold+k1/2;
u2x=reconuxxxp(u2,Nx,dx);
k2=dt*resJSC1(uold,u2x,cc);
u3=uold+k2;
u3x=reconuxxxp(u3,Nx,dx);
k3=dt*resJSC1(uold,u3x,cc);
unew=uold+(k0+2*k1+2*k2+k3)/6;
elseif IRK==2
u0x=reconuxp(uold,Nx,dx);
u0xx=reconuxp(u0x,Nx,dx);
u0xxx=reconuxp(u0xx,Nx,dx);
u2x=reconuxp(uold.^2,Nx,dx);
ko=dt*resJSC1(u2x,u0xxx,cc);
u1=uold+k0;
u1x=reconuxp(u1,Nx,dx);
uxx=reconuxp(u1x,Nx,dx);
uxxx=reconuxp(uxx,Nx,dx);
u2x=reconuxp(u1.^2,Nx,dx);
k1=dt*resJSC1(u2x,uxxx,cc);
u2=u1+k1;
unew=(uold+u2)/2;
end
uold=unew;
eps=0.3*dt;
if abs(k*dt-0.25)<eps|abs(k*dt-0.5)<eps|abs(k*dt-0.75)<eps|abs(k*dt-1)<eps disp(k),
[u0ex]=feval(@JSC1,XX,k*dt,cc);
subplot(2,1,1),plot(XX,unew,'b-',XX,uoex,'g.');
xlabel('x');ylabel('u(x,t)');
hold on
subplot(2,1,2),
if abs(k*dt-0.25)<eps,plot(XX,abs(u0ex-unew),'r-');end
if abs(k*dt-0.5)<eps,plot(XX,abs(u0ex-unew),'r-.');end
if abs(k*dt-0.75)<eps,plot(XX,abs(u0ex-unew),'r--');end
if abs(k*dt-1)<eps,plot(xx,abs(u0ex-unew),'r:');end
xlabel('x');ylabel('inf error');
hold on
end
end
function uu=JSC1(x,ap,cc)
uu=sin(cc*(x+ap));
function uu=resJSC1(u1,u2,cc)
uu=-cc*u2;
function ux=reconuxxxp(u,N,h)
IEX=0;
if IEX==4
a=2;b=-1;
for i=1:N
if i==1
tmp=b*(u(i+3)-3*u(i+1)+3*u(N-1)-u(N-3))/4+a*(u(i+2)-2*u(i+1)+2*u(N-1)-u(N-2));
else if i==2
tmp=b*(u(i+3)-3*u(i+1)+3*u(i-1)-u(N-2))/4+a*(u(i+2)-2*u(i+1)+2*u(i-1)-u(N-1));
else if i==3
tmp=b*(u(i+3)-3*u(i+1)+3*u(i-1)-u(N-1))/4+a*(u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2));
else if i==(N-2)
tmp=b*(u(2)-3*u(i+1)+3*u(i-1)-u(i-3))/4+a*(u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2));
else if i==(N-1)
tmp = b*(u(3)-3*u(i+1)+3*u(i-1)-u(i-3))/4+a*(u(2)-2*u(i+1)+2*u(i-1)-u(i-2));
else if i==N
tmp=b*(u(4)-3*u(2)+3*u(i-1)-u(i-3))/4+a*(u(3)-2*u(2)+2*u(i-1)-u(i-2));
else
tmp=b*(u(i+3)-3*u(i+1)+3*u(i-1)-u(i-3))/4+a*(u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2));
end
ux(i)=tmp/(2*h^3);
end
return
end
if IEX==6
a=-488/240;b=338/240;c=-72/240;d=7/240;
for i=1:N
if i==1
tmp=a*(u(i+1)-u(N-1))+b*(u(i+2)-u(N-2))+c*(u(i+3)-u(N-3))+d*(u(i+4)-u(N-4));
else if i==2
tmp=a*(u(i+1)-u(i-1))+b*(u(i+2)-u(N-1))+c*(u(i+3)-u(N-2))+d*(u(i+4)-u(N-3));
else if i==3
tmp=a*(u(i+1)-u(i-1))+b*(u(i+2)-u(i-2))+c*(u(i+3)-u(N-1))+d(u(i+4)-u(N-2));
else if i==4
tmp=a*(u(i+1)-u(i-1))+b*(u(i+2)-u(i-2))+c*(u(i+3)-u(i-3))+d*(u(i+4)-u(N-1));
else if i==(N-3)
tmp=a*(u(i+1)-u(i-1))+b*(u(i+2)-u(i-2))+c*(u(i+3)-u(i-3))+d(u(2)-u(i-4));
else if i==(N-2)
tmp=a*(u(i+1)-u(i-1))+b*(u(i+2)-u(i-2))+c*(u(2)-u(i-3))+d*(u(3)-u(i-4));
else if i==(N-1)
tmp=a*(u(i+1)-u(i-1))+b*(u(2)-u(i-2))+c*(u(3)-u(i-3))+d*(u(4)-u(i-4));
else if i==N
tmp=a*(u(2)-u(i-1))+b*(u(3)-u(i-2))+c*(u(4)-u(i-3))+d*(u(5)-u(i-4));
else
tmp=a*(u(i+1)-u(i-1))+b*(u(i+2)-u(i-2))+c*(u(i+3)-u(i-3))+d*(u(i+4)-u(i-4));
end
ux(i)=tmp/(h^3);
end
return
end
ISC=6;
if ISC==4
alfa=15/32;aa=2;bb=2*alfa-1;
else if ISC==6
alfa=7/16;aa=2;bb=-1.0/8;
end
amat=zeros(N,N);
B=zeros(N,1);
for i=1:N
if i==1
amat(1,1)=1;amat(1,2)=alfa;amat(1,N-1)=alfa;
B(i)=bb*(u(i+3)-3*u(i+1)+3*u(N-1)-u(N-3))/4+aa*(u(i+2)-2*u(i+1)+2*u(N-1)-u(N-2));
else if i==2
amat(2,1)=alfa;amat(2,2)=1;amat(2,3)=alfa;
B(i)=bb*(u(i+3)-3*u(i+1)+3*u(i-1)-u(N-2))/4+aa*(u(i+2)-2*u(i+1)+2*u(i-1)-u(N-1));
else if i==3
amat(i,i-1)=alfa;amat(i,i)=1;amat(i,i+1)=alfa;
B(i)=bb*(u(i+3)-3*u(i+1)+3*u(i-1)-u(N-1))/4+aa*(u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2));
else if i==N-2
amat(i,i-1)=alfa;amat(i,i)=1;amat(i,i+1)=alfa;
B(i)=bb*(u(2)-3*u(i+1)+3*u(i-1)-u(i-3))/4+aa*(u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2));
else if i==N-1
amat(i,i-1)=alfa;amat(i,i)=1;amat(i,i+1)=alfa;
B(i)=bb*(u(3)-3*u(i+1)+3*u(i-1)-u(i-3))/4+aa*(u(2)-2*u(i+1)+2*u(i-1)-u(i-2));
else if i==N
amat(N,1)=-1;amat(N,N)=1;
B(i)=0
else
amat(i,i-1)=alfa;amat(i,i)=1;amat(i,i+1)=alfa;
B(i)=bb*(u(i+3)-3*u(i+1)+3*u(i-1)-u(i-3))/4+aa*(u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2));
end
B(i)=B(i)/(2*h^3);
ux=(amat\B)';
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
  1 comentario
KSSV
KSSV el 15 de Sept. de 2021
What error you are getting? Your ocde has got lot of loops. I think there is a lot of room to make your code effective.

Iniciar sesión para comentar.

Iniciar sesión para responder a esta pregunta.

Respuestas (1)

Walter Roberson
Walter Roberson el 15 de Sept. de 2021
The function you are missing is on PDF page 108 of https://www.researchgate.net/file.PostFileLoader.html?id=5696da346143257e508b4567&assetKey=AS%3A317578013020160%401452727906601
  2 comentarios
Star Strider
Star Strider el 15 de Sept. de 2021
+1
... because I have absolutely no idea how you figured that out!
.
Walter Roberson
Walter Roberson el 15 de Sept. de 2021
just googled reconuxp and looked at the first few results.

Iniciar sesión para comentar.

Categorías

Más información sobre Mathematics en Help Center y File Exchange.

Etiquetas

Productos

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by