solve解方程遇到问题（已看过别人的帖子没有思路，所以求助）

ziding1763

首先说明：已看过论坛中有关solve的帖子，问题没有解决，所以请大家帮忙。

问题描述：以矩形的一个顶点为轴，旋转一个矩形，旋转角度为phi.由三个小程序共同实现。
主程序
function demo
border_x0=3;
border_x1=4;%固定矩形
border_y0=6;
border_y1=9;
h0=0.4;h=0.05;phi=pi/6;
fix1=mygetpoint(border_x0,border_x1,border_y0,border_y1,phi,h0,h)

调用程序
function fix=mygetpoint(border_x0,border_x1,border_y0,border_y1,phi,h0,h)
% This function is used to find the fixed point in the delaunay mesh of the
% deleted boundary. border_x0,border_x1,border_y0,border_y1 are the
% displacement of the area. h0 is the size of edge length. h is the length
% of the x-direction step.
x=border_x0:h:border_x1;
y=border_y0:h:border_y1;
x1=border_x1:-0.05:border_x0;
y1=border_y1:-0.4:border_y0;
centre_x=0.5*(border_x0+border_x1);
centre_y=0.5*(border_y0+border_y1);
p1=[x(:) ones(length(x),1)*border_y0];%按照每个边，固定出矩形的一些点
p2=[ones(length(y),1)*border_x1  y(:)];
p3=[x1(:) ones(length(x1),1)*border_y1];
p4=[ones(length(y1),1)*border_x0 y1(:)];

p=[p1;p2;p3;p4];
axis([0 20 0 16])
plot(p(:,1),p(:,2))
p=myprotate0(p,phi);
hold on
p=[p(:,1) p(:,2)];
plot(p(:,1),p(:,2))
fix=p;

另一个调用程序
function p=myprotate0(p,phi)
[M,N]=size(p);
for i=2:M
    d0=sqrt((p(i,1)-p(1,1))^2+(p(i,2)-p(1,2))^2);%矩形上的点到固定定点的距离
    [x,y]=solve('2*d0^2-sqrt((x-p(i,1))^2+(y-p(i,2))^2)=2*d0^2*cos(phi)','(x-p(1,1))^2+(y-p(1,2))^2=d0^2','x,y')%解出旋转后点的坐标。
    x=subs(x);y=subs(y);%这是为了solve 中的值为sym型的，而想把它转为数值型，问题出现的地方。
    p(i,1)=x;p(i,2)=y;   
end


出现的错误：
x =

[ -1/2*(-1/(-8*p(i,2)*p(1,2)+4*p(1,2)^2+4*p(i,1)^2+4*p(1,1)^2-8*p(i,1)*p(1,1)+4*p(i,2)^2)*(16*d0^4*p(1,2)+4*p(1,2)^3-4*p(1,2)*d0^2+4*p(i,2)^3+32*p(i,2)*d0^4*cos(phi)+4*p(i,1)^2*p(i,2)+4*p(i,1)^2*p(1,2)-16*p(i,2)*d0^4+4*p(i,2)*p(1,1)^2-8*p(1,1)*p(i,1)*p(i,2)-4*p(i,2)*p(1,2)^2+4*p(i,2)*d0^2-4*p(i,2)^2*p(1,2)+4*p(1,1)^2*p(1,2)-8*p(1,1)*p(i,1)*p(1,2)-32*d0^4*cos(phi)*p(1,2)-16*p(i,2)*d0^4*cos(phi)^2+16*d0^4*cos(phi)^2*p(1,2)+4*(8*p(1,2)*d0^2*p(1,1)*p(i,1)*p(i,2)+8*d0^4*p(1,2)^2*p(i,1)^2+8*d0^4*p(1,2)^2*p(1(1,1)*p(i,1)*p(i,2)*d0^4*cos(phi)^2*p(1,2)-16*p(i,1)^2*p(i,2)*d0^4*cos(phi)^2*p(1,2)+32*p(i,2)^2*d0^4*cos(phi)*p(1,1)*p(i,1)+32*p(i,2)*d0^4*cos(phi)*p(1,1)^2*p(1i,2)-4*p(1,2)*d0^2*p(i,2)*p(1,1)^2+6*p(1,1)^5*p(i,1)-16*p(1,1)^2*d0^8+8*p(1,1)^2*d0^6+2*p(1,1)^4*d0^2+8*p(1,1)^4*d0^4-2*p(1,1)^4*p(1,2)^2-p(1,2)^4*p(1,1)^2-p(i,2)^4*p(i,1)^2-p(i,2……………………………………
y =

[ 1/2/(-8*p(i,2)*p(1,2)+4*p(1,2)^2+4*p(i,1)^2+4*p(1,1)^2-8*p(i,1)*p(1,1)+4*p(i,2………
Conversion to double from sym is not possible.
Error in ==> d:\MATLAB6p5\work\myprotate0.m
On line 8  ==>     p(i,1)=x;p(i,2)=y;   
Error in ==> d:\MATLAB6p5\work\mygetpoint.m
On line 22  ==> p=myprotate0(p,phi);
Error in ==> d:\MATLAB6p5\work\demo.m
On line 7  ==> fix1=mygetpoint(border_x0,border_x1,border_y0,border_y1,phi,h0,h)

上面的x,y还有很长的公式，这里只是说明一下，可以看出来在其中计算的本来是调用时可以确定的数值p(i,1),p(i,2),p(1,1),p(1,2),d0等这些都是还用原有的符号的行使给出。不能带到原来的式子中计算数值解吗？

请大家帮忙改程序错误，或者有更好的实现方法，请提供。谢谢。

[ 本帖最后由 eight 于 2007-8-12 12:11 编辑 ]