I am having error in this programme. i gave the call function also.
Mostrar comentarios más antiguos
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
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.
Respuestas (1)
Walter Roberson
el 15 de Sept. de 2021
1 voto
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
el 15 de Sept. de 2021
+1
... because I have absolutely no idea how you figured that out!
.
Walter Roberson
el 15 de Sept. de 2021
just googled reconuxp and looked at the first few results.
Categorías
Más información sobre Programming en Centro de ayuda y File Exchange.
Productos
el 15 de Sept. de 2021
el 15 de Sept. de 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
