坐标变换程序运行得不到结果

kerb · 发表于 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)

你检查一下是否搞混了

hipeilei · 发表于 2015-11-13 23:10:24

本帖最后由 hipeilei 于 2015-11-14 00:12 编辑

kerb 发表于 2015-11-13 00:16
另外用迭代法求解非线性方程组应该是这样

多谢指点，今天考了一天的试，到现在才看到，感谢关注。我感觉也是雅可比矩阵这里出了问题，可就是不知道到底错在什么地方了，我用maple检查了一下我的雅可比矩阵，和我代码里面的是一样的，感觉没求错。下面是maple求导的结果：

这是代码里面的式子，其中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；那就说明是我的雅可比矩阵求的有问题

答主应该也是力学的吧，方便qq交流吗？我的qq：913111805

hipeilei · 发表于 2015-11-13 23:37:45

kerb 发表于 2015-11-12 23:42
好像可以解出两个坐标系变换，稍微麻烦一点

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

pasuka · 发表于 2015-11-15 16:50:46

hipeilei 发表于 2015-11-13 23:37
这样解出来可能会有两个根，取舍问题不好判断，有可能两个点都在规则单元里。 ...

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

kerb · 发表于 2015-11-18 15:50:07

你手头有maple，直接解析求出(\xi,\eta)的表达式，然后稍加判断就行了

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

水王勋章

元老勋章

热心勋章