去網(wǎng)上抄了個(gè)畫(huà)圖的函數(shù),但是報(bào)錯(cuò)了,求各位大神解答#sol_equ用來(lái)解六次方程sol_equ<-function(coe1,coe2,coe3,coe4,coe5,coe6,coe7,b1,b2,e){#定義了一個(gè)六次函數(shù)cubiccubic<-function(x){return(coe1*x^6+coe2*x^5+coe3*x^4+coe4*x^3+coe5*x^2+coe6*x+coe7)}k=0while(2>1){if(k>15){#設(shè)置最大循環(huán)次數(shù)是15return(NA)}k=k+1y1<-cubic(b1)y2<-cubic(b2)b0<-(b1+b2)/2if(cubic(b0)==0){return(b0)}if(y1==0){return(b1)}if(y2==0){return(b2)}y0<-cubic(b0)if(y1*y0<0){b2<-b0}if(y2*y0<0){b1<-b0}if(abs(b1-b2)<e){result<-((b1+b2)/2)break}}return(result)}#設(shè)置y,z 坐標(biāo)值y<-seq(-1.5,1.5,0.01)z<-seq(-1.5,1.5,0.01)x<-NULL#下面是對(duì)于每個(gè)(y,z)(共301×301)求得x(x>0)for(i in 1:301){m<-NULLfor(j in 1:301){a<-z[j]^2+9*y^2/4-1b<–z[j]^3c<–9*y^2*z[j]^3/80r<-sol_equ(1,0,3*a,0,3*a^2+b,0,a^3+c,0,2,0.01)m<-c(m,r)}x<-cbind(x,m)}xlibrary(rgl)#當(dāng)x>0時(shí)作圖persp3d(y,z,t(x),theta=60,phi=30,col=’red’)#當(dāng)x<0時(shí)作圖persp3d(y,z,-t(x),theta=60,phi=30,col=’red’,add=T)報(bào)錯(cuò)如圖
R語(yǔ)言函數(shù)問(wèn)題,求解
對(duì)方正在青青草原抓羊
2016-10-17 21:46:17