[已解决] MATLAB新手,运行出现索引超出数组元素的数目,能帮忙看一下吗?

[复制链接]
wyzwyz123 发表于 2021-10-11 09:20:24



  1. sol=pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t)
  2. m=2;
  3. x = linspace(0,1,50);
  4. t = linspace(0.1,pi,50);
  5. [pL,qL,pR,qR]=feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),y0,t(1),varargin{:});

  6. [pL,qL,pR,qR]=feval(bc,x(1),u(:,1),x(nx),u(:,nx),u,tnow,varargin{:});

  7. surf(x,t,sol(:,:,1))

  8. figure;

  9. kdf=[d0,kf];
  10. function [c,f,s]=pdex4pde(x,t,u,dudx)
  11. global kdf
  12. c=1/(kdf(1))/exp(1*u./(0.737*0.2079^(1/1.41)));
  13. f=dudx;
  14. s=(1/(0.737*0.2079^(1/1.41)))*(dudx.^2);
  15. end

  16. function u0=pdex4ic(x)
  17. u0=0;
  18. end


  19. function [pl,ql,pr,qr]=pdex4bc(xl,ul,xr,ur,u,t)


  20. global kdf
  21. u=abs(u);
  22. pl=0;
  23. ql=0;
  24. n=1000;
  25. r=0.000593/2;
  26.                     
  27. h=r/(n-1);
  28. s=zeros (1000,1);
  29. a=0;
  30. for i=1:1:n
  31.     if i==1
  32.       s(i,1)=4*pi/3*(h/2)^3.*u(:,i);
  33.       a=a+s(i,1);
  34.     end
  35.     if i==n
  36.       s(i,1)=(4*pi/3*r^3-4*pi/3*(r-h/2)^3).*u(:,i);
  37.       a=a+s(i,1);
  38.     end
  39.     if i>1 & i<n   
  40.      s(i,1)=(4*pi/3*(h*(i-1)+h/2)^3-4*pi/3*(h*(i-1)-h/2)^3).* u(:,i);
  41.      a=a+s(i,1);
  42.     end
  43. end
  44. ss=(100*0.2079-30*a/(4/3*pi*(r)^3))/(100-30*a/(4/3*pi*(r)^3))-(u(:,n)/0.737)^1.41;  
  45. pr=kdf(2)*ss;
  46. qr=-kdf(1)*(exp(1*ul/(0.737*0.2079^(1/1.41))));
  47. end

复制代码


最佳答案


20141303 发表于 2021-10-11 09:46:11
仅供参考,首先该代码中sol=pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t)这句代码应在t = linspace(0.1,pi,50);之后,其次在子函数中global kdf声明了全局变量,但在主函数中也应有global kdf声明全局变量
回复此楼

6 条回复


20141303 发表于 2021-10-11 09:46:11
仅供参考,首先该代码中sol=pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t)这句代码应在t = linspace(0.1,pi,50);之后,其次在子函数中global kdf声明了全局变量,但在主函数中也应有global kdf声明全局变量
回复此楼

wyzwyz123 发表于 2021-10-12 09:37:14
20141303 发表于 2021-10-11 09:46
仅供参考,首先该代码中sol=pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t)这句代码应在t = linspace(0.1,pi,50 ...

您好,我对kdf进行定义了,为什么还是不能运行?

  1. function main
  2. global kdf
  3. d0 = 1;
  4. kf = 1;
  5. kdf=[d0,kf];
  6. m=2;
  7. x = linspace(0,1,50);
  8. t = linspace(0.1,pi,50);
  9. sol=pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
  10. [pL,qL,pR,qR]=feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),y0,t(1));
  11. [pL,qL,pR,qR]=feval(bc,x(1),u(:,1),x(nx),u(:,nx),u,tnow);
  12. figure;
  13. surf(x,t,sol(:,:,1))
  14. end
  15. function [c,f,s]=pdex4pde(x,t,u,dudx)
  16. global kdf
  17. c=1/(kdf(1))/exp(1*u./(0.737*0.2079^(1/1.41)));
  18. f=dudx;
  19. s=(1/(0.737*0.2079^(1/1.41)))*(dudx.^2);
  20. end
  21. function u0=pdex4ic(x)
  22. u0=0;
  23. end

  24. function [pl,ql,pr,qr]=pdex4bc(xl,ul,xr,ur,u,t)
  25. global kdf
  26. u=abs(u);
  27. pl=0;
  28. ql=0;
  29. n=1000;
  30. r=0.000593/2;
  31. h=r/(n-1);
  32. s=zeros (1000,1);
  33. a=0;
  34. for i=1:1:n
  35.     if i==1
  36.       s(i,1)=4*pi/3*(h/2)^3.*u(:,i);
  37.       a=a+s(i,1);
  38.     end
  39.     if i==n
  40.       s(i,1)=(4*pi/3*r^3-4*pi/3*(r-h/2)^3).*u(:,i);
  41.       a=a+s(i,1);
  42.     end
  43.     if i>1 && i<n   
  44.      s(i,1)=(4*pi/3*(h*(i-1)+h/2)^3-4*pi/3*(h*(i-1)-h/2)^3).* u(:,i);
  45.      a=a+s(i,1);
  46.     end
  47. end
  48. ss=(100*0.2079-30*a/(4/3*pi*(r)^3))/(100-30*a/(4/3*pi*(r)^3))-(u(:,n)/0.737)^1.41;  
  49. pr=kdf(2)*ss;
  50. qr=-kdf(1)*(exp(1*ul/(0.737*0.2079^(1/1.41))));
  51. end
复制代码

20141303 发表于 2021-10-12 09:45:09
仅供参考,子函数pdex4bc输入参数u为单个数字,故索引u(:,i)等报错,将其等改为u试试

wyzwyz123 发表于 2021-10-12 09:51:29
20141303 发表于 2021-10-12 09:45
仅供参考,子函数pdex4bc输入参数u为单个数字,故索引u(:,i)等报错,将其等改为u试试
...

是把abs()去掉吗?

20141303 发表于 2021-10-12 09:52:55
仅供参考
  1. % function main
  2. global kdf
  3. d0 = 1;
  4. kf = 1;
  5. kdf=[d0,kf];
  6. m=2;
  7. x = linspace(0,1,50);
  8. t = linspace(0.1,pi,50);
  9. sol=pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
  10. [pL,qL,pR,qR]=feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),y0,t(1));
  11. [pL,qL,pR,qR]=feval(bc,x(1),u(:,1),x(nx),u(:,nx),u,tnow);
  12. figure;
  13. surf(x,t,sol(:,:,1))
  14. % end
  15. function [c,f,s]=pdex4pde(x,t,u,dudx)
  16. global kdf
  17. c=1/(kdf(1))/exp(1*u./(0.737*0.2079^(1/1.41)));
  18. f=dudx;
  19. s=(1/(0.737*0.2079^(1/1.41)))*(dudx.^2);
  20. end
  21. function u0=pdex4ic(x)
  22. u0=0;
  23. end

  24. function [pl,ql,pr,qr]=pdex4bc(xl,ul,xr,ur,u,t)
  25. global kdf
  26. u=abs(u);
  27. pl=0;
  28. ql=0;
  29. n=1000;
  30. r=0.000593/2;
  31. h=r/(n-1);
  32. s=zeros (1000,1);
  33. a=0;
  34. for i=1:1:n
  35.     if i==1
  36.       s(i,1)=4*pi/3*(h/2)^3.*u;
  37.       a=a+s(i,1);
  38.     end
  39.     if i==n
  40.       s(i,1)=(4*pi/3*r^3-4*pi/3*(r-h/2)^3).*u;
  41.       a=a+s(i,1);
  42.     end
  43.     if i>1 && i<n   
  44.      s(i,1)=(4*pi/3*(h*(i-1)+h/2)^3-4*pi/3*(h*(i-1)-h/2)^3)*u;
  45.      a=a+s(i,1);
  46.     end
  47. end
  48. ss=(100*0.2079-30*a/(4/3*pi*(r)^3))/(100-30*a/(4/3*pi*(r)^3))-(u/0.737)^1.41;  
  49. pr=kdf(2)*ss;
  50. qr=-kdf(1)*(exp(1*ul/(0.737*0.2079^(1/1.41))));
  51. end
复制代码

wyzwyz123 发表于 2021-10-12 10:31:13

您好,我把
  1. function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
复制代码

通过
  1. [pL,qL,pR,qR]=feval(bc,x(1),u(:,1),x(nx),u(:,nx),u,tnow);
复制代码

改为
  1. function [pl,ql,pr,qr]=pdex4bc(xl,ul,xr,ur,u,t)
复制代码
可行吗?它总显示Untitled4>pdex4bc ,Untitled4是我整个代码的命名
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表