Fortran Coder

楼主: hipeilei

[空间几何] 坐标变换程序运行得不到结果

[复制链接]

59

帖子

2

主题

0

精华

大师

F 币
810 元
贡献
476 点
发表于 2015-11-13 00:16:35 | 显示全部楼层
本帖最后由 kerb 于 2015-11-13 00:22 编辑

另外用迭代法求解非线性方程组应该是这样
F(\xi,\eta)=\left\{\begin{array}{c}f_1(\xi,\eta)=a_x+b_x\xi+c_x\eta+d_x\xi\eta-x\\f_2(\xi,\eta)=a_y+b_y\xi+c_y\eta+d_y\xi\eta-y\end{array}\right.

F(0+\delta)=F(0)+\nabla F\cdot\delta+O(\delta^2),\delta=(\Delta\xi,\Delta\eta)所以那个雅克比矩阵是对\xi,\eta求导

也就是:
J=\left(\begin{array}{cc}b_x+d_x\eta & c_x+d_x\xi\\b_y+d_y\eta & c_y+d_y\xi\end{array}\right)

你检查一下是否搞混了

18

帖子

3

主题

0

精华

入门

F 币
85 元
贡献
48 点
 楼主| 发表于 2015-11-13 23:10:24 | 显示全部楼层
本帖最后由 hipeilei 于 2015-11-14 00:12 编辑
kerb 发表于 2015-11-13 00:16
另外用迭代法求解非线性方程组应该是这样

多谢指点,今天考了一天的试,到现在才看到,感谢关注。我感觉也是雅可比矩阵这里出了问题,可就是不知道到底错在什么地方了,我用maple检查了一下我的雅可比矩阵,和我代码里面的是一样的,感觉没求错。下面是maple求导的结果: 35435.png
这是代码里面的式子,其中m相当于kesi,n相当于yita,a表示x对kesi求偏导,b表示y对kesi求偏导,c表示x对yita求偏导,d表示y对yita求偏导:
subroutine Jacobi(x1,y1,x2,y2,x3,y3,x4,y4,m,n,a,b,c,d)
implicit none
real:: a,b,c,d,m,n
real,dimension(2,2)::J
real:: x1,y1,x2,y2,x3,y3,x4,y4
a=(x1*(-1+n)+x2*(1-n)+x3*(1+n)+x4*(-1-n))/4
b=(y1*(-1+n)+y2*(1-n)+y3*(1+n)+y4*(-1-n))/4
c=(x1*(-1+m)+x2*(-1-m)+x3*(1+m)+x4*(1-m))/4
d=(y1*(-1+m)+y2*(-1-m)+y3*(1+m)+y4*(1-m))/4

  J(1,1)=d/(a*d-b*c)
  J(1,2)=-b/(a*d-b*c)
  J(2,1)=-c/(a*d-b*c)
  J(2,2)=a/(a*d-c*b)
write(*,*),"***********This is Jacobi matrix!*******************"
write(*,"(4f10.2)")J(1,1),J(1,2),J(2,1),J(2,2)
end subroutine Jacobi
调试的时候发现了一个新问题,就是雅可比矩阵求出来居然是0;那就说明是我的雅可比矩阵求的有问题 Capture.PNG
答主应该也是力学的吧,方便qq交流吗?我的qq:913111805



18

帖子

3

主题

0

精华

入门

F 币
85 元
贡献
48 点
 楼主| 发表于 2015-11-13 23:37:45 | 显示全部楼层
kerb 发表于 2015-11-12 23:42
好像可以解出两个坐标系变换,稍微麻烦一点

这样解出来可能会有两个根,取舍问题不好判断,有可能两个点都在规则单元里。

490

帖子

4

主题

0

精华

大宗师

F 币
3298 元
贡献
1948 点

水王勋章元老勋章热心勋章

发表于 2015-11-15 16:50:46 | 显示全部楼层
hipeilei 发表于 2015-11-13 23:37
这样解出来可能会有两个根,取舍问题不好判断,有可能两个点都在规则单元里。 ...

碰上二义性吧?一般四边形处理云图或等值线会遇到,解决办法有很多
最直观的办法就是:四边形加条对角线切成2个三角形,假定点q不在四条边上,个么q要么在两个三角形内部,要么在对角线上,等参变换后亦如此
至于如何判断q是否在三角形内部,有个帖子已经介绍过了

59

帖子

2

主题

0

精华

大师

F 币
810 元
贡献
476 点
发表于 2015-11-18 15:50:07 | 显示全部楼层
你手头有maple,直接解析求出(\xi,\eta)的表达式,然后稍加判断就行了
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-3-29 19:42

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表