| 查看: 2908 | 回复: 2 | ||
[求助]
求助一个matlab程序的问题 已有1人参与
|
|
m文件如下 function R = R4(f,g,m,n,a,b,xa,ya,ua,va,N) f=@(u)(u); g=@(v)(v); m=@(x,y,u)(5.487*(10^20)*(16*x*y/pi+(1-4*y)*(1/pi+log(4*x)))*x/(sqrt(x*x+y*y)-3468*u)); n=@(x,y,v)(5.487*(10^20)*(16*x*y/pi+(1-4*y)*(1/pi+log(4*x)))*y/(sqrt(x*x+y*y)-3468*(v-1))); h=(b-a)/N; T=zeros(1,N+1); X=zeros(1,N+1); Y=zeros(1,N+1); U=zeros(1,N+1); V=zeros(1,N+1); for n=1:N K1=feval(f,U(n)); K2=feval(f,U(n)+h/2*K1); K3=feval(f,U(n)+h/2*K2); K4=feval(f,U(n)+h*K3); X(n+1)=X(n)+h*(K1+2*K2+2*K3+K4)/6; L1=feval(g,V(n)); L2=feval(g,V(n)+h/2*L1); L3=feval(g,V(n)+h/2*L2); L4=feval(g,V(n)+h*L3); Y(n+1)=Y(n)+h*(L1+2*L2+2*L3+L4)/6; M1=feval(m,X(n),Y(n),U(n)); N1=feval(n,X(n),Y(n),V(n)); M2=feval(m,X(n)+h/2*K1,Y(n)+h/2*L1,U(n)+h/2*M1); N2=feval(n,X(n)+h/2*K1,Y(n)+h/2*L1,V(n)+h/2*N1); M3=feval(m,X(n)+h/2*K2,Y(n)+h/2*L2,U(n)+h/2*M2); N3=feval(n,X(n)+h/2*K2,Y(n)+h/2*L2,V(n)+h/2*N2); M4=feval(m,X(n)+h*K3,Y(n)+h*L3,U(n)+h*M3); N4=feval(n,X(n)+h*k3,Y(n)+h*L3,V(n)+h*N3); U(n+1)=U(n)+h*(M1+2*M2+2*M3+M4); V(n+1)=V(n)+h*(N1+2*N2+2*N3+N4); R=[X' Y' U' V']; end 然后命令行输入 R4('f','g','m','n',0,250,0,0,0,1,250)后出错 错误使用 feval 参数必须包含字符矢量或函数句柄。 出错 R4 (line 24) N1=feval(n,X(n),Y(n),V(n)); |
回复此楼» 猜你喜欢
全日制（定向）博士
已经有5人回复
-
假如你的研究生提出不合理要求
已经有10人回复
-
萌生出自己或许不适合搞科研的想法，现在跑or等等看？
已经有4人回复
-
Materials Today Chemistry审稿周期
已经有4人回复
-
参与限项
已经有3人回复
-
实验室接单子
已经有4人回复
-
对氯苯硼酸纯化
已经有3人回复
-
求助:我三月中下旬出站，青基依托单位怎么办？
已经有12人回复
-
所感
已经有4人回复
-
要不要辞职读博？
已经有7人回复
2楼2017-05-08 12:21:38
1314168apple
金虫 (知名作家)
- 应助: 68 (初中生)
- 金币: 677
- 红花: 12
- 帖子: 6872
- 在线: 1462.3小时
- 虫号: 287760
- 注册: 2006-10-21
- 专业: 色谱分析
【答案】应助回帖
感谢参与，应助指数 +1
|
很多错误， 1、函数如何定义都不知道。f、g、m、n属于函数，不能够在输入处。 2、定义n为函数，在循环的地方再次使用n使之变为数值。 3、xa，ya，ua，va不知道用在那里？ 4、大小写不注意 K3、k3。 5、T也不知所谓。 6、f、g的函数也多此一举(我这里不修改了)。 X(1)、Y(1)、U(1)、V(1)的值的为零，这不符合。 看你的命令好像是希望计算微分方程，用的是4-5阶Runge-Kutta 。其实可以用自带命令ode45. function R = R4(a,b,N) f=@(u)(u); g=@(v)(v); m=@(x,y,u)(5.487*(10^20)*(16*x*y/pi+(1-4*y)*(1/pi+log(4*x)))*x/(sqrt(x*x+y*y)-3468*u)); k=@(x,y,v)(5.487*(10^20)*(16*x*y/pi+(1-4*y)*(1/pi+log(4*x)))*y/(sqrt(x*x+y*y)-3468*(v-1))); h=(b-a)/N; T=zeros(1,N+1); X=zeros(1,N+1); Y=zeros(1,N+1); U=zeros(1,N+1); V=zeros(1,N+1); for n=1:N K1=feval(f,U(n)); K2=feval(f,U(n)+h/2*K1); K3=feval(f,U(n)+h/2*K2); K4=feval(f,U(n)+h*K3); X(n+1)=X(n)+h*(K1+2*K2+2*K3+K4)/6; L1=feval(g,V(n)); L2=feval(g,V(n)+h/2*L1); L3=feval(g,V(n)+h/2*L2); L4=feval(g,V(n)+h*L3); Y(n+1)=Y(n)+h*(L1+2*L2+2*L3+L4)/6; M1=feval(m,X(n),Y(n),U(n)); N1=feval(k,X(n),Y(n),V(n)); M2=feval(m,X(n)+h/2*K1,Y(n)+h/2*L1,U(n)+h/2*M1); N2=feval(k,X(n)+h/2*K1,Y(n)+h/2*L1,V(n)+h/2*N1); M3=feval(m,X(n)+h/2*K2,Y(n)+h/2*L2,U(n)+h/2*M2); N3=feval(k,X(n)+h/2*K2,Y(n)+h/2*L2,V(n)+h/2*N2); M4=feval(m,X(n)+h*K3,Y(n)+h*L3,U(n)+h*M3); N4=feval(k,X(n)+h*K3,Y(n)+h*L3,V(n)+h*N3); U(n+1)=U(n)+h*(M1+2*M2+2*M3+M4); V(n+1)=V(n)+h*(N1+2*N2+2*N3+N4); end R=[X' Y' U' V']; R4(0,250,250) ans = 0 0 0 0 0 0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 。。。。。。。。。。。。。 |
3楼2017-05-11 09:13:55