-1

这是我想扩展和派生的函数的乳胶代码：

$f(x) = \sum\limits_{i=1}^n \left(x_i + \frac{h}{2}(1-ih)\sum\limits_{j=1}^i jh(x_j + jh +1)^3 + \frac{h}{2}ih\sum\limits_{j=i+1}^n (1-jh)(x_j+jh+1)^3 \right)^2$

在此处输入图像描述

我计划n=6通过deriv(expression(myfunction),c('x1', 'x2','x3','x4','x5','x6')). 我不知道如何在不扩展函数的情况下计算导数，但如果有办法请告诉我。

我的问题始于尝试仅在一个字符串中扩展函数时，因为部分表达式会重复，因此将不胜感激：

n<-6
myexpr <- sapply(1:n, function(i) paste( paste('x',i,sep=''),paste('+h/2*(1-',i,'*h)*(',sep=''),
                           sapply(1:i,function(j) paste('(j*h*(',paste('x',j,sep=''),'j*h+1)^3)',collapse='+')),
                           paste('+h/2*',i,'*h*(',sep=''),
                           sapply((i+1):n,function(j) paste('((1-j*h)*(',paste('x',j,sep=''),'+j*h+1)^3)',collapse='+'))
                           ,collapse='+'))


deriv(expression(myexpr),c('x1', 'x2','x3','x4','x5','x6'))
4

2 回答 2

2

我认为您应该尝试一些更简单的方法来处理表达式。在您的情况下，只需查看第一个参数，我就可以在您构造这个（大得离谱的）表达式中发现一个错误。

> do.call(deriv, list(expr=parse(text=myexpr[1]), namevec=c('x1') ) )
Error in parse(text = myexpr[1]) : <text>:1:31: unexpected symbol
1: x1 +h*0.5*(1-1*h)*( (j*h*( x1 j
                                  ^
> substr(myexpr[1],1,40)
[1] "x1 +h*0.5*(1-1*h)*( (j*h*( x1 j*h+1)^3) "

因此，您在第二项的扩展中缺少“+”号。

然而，从战略上讲，我会认为 Mathematica 或 Maxima 会是用于此目的的更好平台。

于 2013-10-18T16:24:01.793 回答
0

我被否决了一些我认为可以做到的事情。事实上，这不仅可以做到，而且 R 解决得非常快：

z <- paste(
  sapply(1:6, function(i,n=6) {

    require(MASS)
    h <- fractions(1/(n+1))

    a <- paste('x',i,sep='')
    b <- paste('+.5*',h,'*(1-',i,'*',h,')*(',sep='')
    cc <- paste(sapply(1:i,function(j) paste(j,'*',h,'*(',paste('x',j,sep=''),'+',j,'*',h,'+1)^3',sep='')),collapse='+')
    d <- paste(a,b,cc,')')

    e <- if (i < 6) paste('+.5*',h,'*(',i,'*',h,')*(',sep='') else ''
    f <- if (i < 6) paste(sapply((i+1):n,function(j) paste('(1-',j,'*',h,')*(',paste('x',j,sep=''),'+',j,'*',h,'+1)^3',sep='')),collapse='+') else ''
    g <- paste(e,f,ifelse(i!=6,')',""))

    paste('(',d,g,')^2',sep='')

  }
  ),
  collapse='+')

deriv(parse(text=z),c('x1','x2','x3','x4','x5','x6'))
于 2013-10-23T12:13:27.393 回答