[已解决] 符号矩阵积分转成数值积分

[复制链接]
congyt 发表于 2022-4-12 09:00:11
一个矩阵,元素都是符号表达式,因为过于复杂导致计算速度太慢,所以我想转成数值积分,但是出现了错误,请大家指点迷津,谢谢!


代码如下:

>> syms x y

A=[x+1 x-1; y+2*x y-1];
B=[2*x+1 3*x-1; 2*y+2*x 3*y-1];
C=A*B;

g = matlabFunction(C);

quad2d(@(x,y)g(x,y),0,1,0,1);
错误使用 reshape
元素数不能更改。请使用 [] 作为大小输入之一,以自动计算该维度的适当大小。
出错
symengine>@(x,y)reshape([(x.*2.0+y.*2.0).*(x-1.0)+(x.*2.0+1.0).*(x+1.0),(x.*2.0+1.0).*(x.*2.0+y)+(x.*2.0+y.*2.0).*(y-1.0),(x.*3.0-1.0).*(x+1.0)+(y.*3.0-1.0).*(x-1.0),(x.*3.0-1.0).*(x.*2.0+y)+(y.*3.0-1.0).*(y-1.0)],[2,2])
出错 @(x,y)g(x,y)

出错 quad2d/tensor (第 343 行)
        Z = FUN(X,Y);  NFE = NFE + 1;
出错 quad2d (第 167 行)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);


最佳答案


20141303 发表于 2022-4-12 09:41:39
仅供参考
  1. clear
  2. clc
  3. syms x y

  4. A=[x+1 x-1; y+2*x y-1];
  5. B=[2*x+1 3*x-1; 2*y+2*x 3*y-1];
  6. C=A*B;
  7. for i=1:2
  8.     for j=1:2
  9.         g = eval(['@(x,y)',vectorize(C(i,j))]);
  10.         s(i,j)=integral2(g,0,1,0,1);
  11.     end
  12. end
复制代码
回复此楼

5 条回复


20141303 发表于 2022-4-12 09:41:39
仅供参考
  1. clear
  2. clc
  3. syms x y

  4. A=[x+1 x-1; y+2*x y-1];
  5. B=[2*x+1 3*x-1; 2*y+2*x 3*y-1];
  6. C=A*B;
  7. for i=1:2
  8.     for j=1:2
  9.         g = eval(['@(x,y)',vectorize(C(i,j))]);
  10.         s(i,j)=integral2(g,0,1,0,1);
  11.     end
  12. end
复制代码
回复此楼

congyt 发表于 2022-4-12 09:50:53

谢谢!因为矩阵比较大,用for循环的话还是比较慢,请问如果不用for循环的话怎么做呢?

20141303 发表于 2022-4-12 10:23:13
抱歉,我也不想用for循环,但目前我也没啥好办法

TouAkira 发表于 2022-4-12 11:05:17
congyt 发表于 2022-4-11 21:50
谢谢!因为矩阵比较大,用for循环的话还是比较慢,请问如果不用for循环的话怎么做呢? ...

仍然是有方法的
  1. LB_x = 0; UB_x = 1; LB_y = 0; UB_y = 1;
  2. C_Func = matlabFunction( C, 'vars', [ x, y ] );
  3. s_vec = integral( @( y ) integral( @( x ) C_Func( x, y ), LB_x, UB_x, 'ArrayValued', 1 ), LB_y, UB_y, 'ArrayValued', 1 );
  4. disp( s_vec )
复制代码

可得
          2.33333333333333                      0.75
                       2.5                      1.25
与循环计算所得结果对比,绝对误差在 e-13 数量级,可以忽略。

congyt 发表于 2022-4-19 09:12:24
TouAkira 发表于 2022-4-12 11:05
仍然是有方法的

可得

谢谢!这几天我用的上边那个方法,我会再用您的方法试一下,非常感谢!
您需要登录后才可以回帖 登录 | 注册

本版积分规则

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