应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » MATLAB/Mathematica » 请教matlab用数组给参数矩阵赋值的问题
121/2返回列表12下一页
查看: 2402  |  回复: 11
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

加油大楠楠

铜虫 (小有名气)

[求助] 请教matlab用数组给参数矩阵赋值的问题

syms x y z d1 d2 d3 r l;n=1;%定义矩阵参数
a=1/l*[z+r-d1 0 0;0 -0.866*y-0.5*z+r-d2 0;0 0 0.866*y-0.5*z+r-d3];
b=[x y z+r-d1;x y-0.866*r+0.866*d2 z-0.5*r+0.5*d2;x y+0.866*r-0.866*d3 z-0.5*r+0.5*d3];
c=inv(a);
d=c*b;%d为我需要的比矩阵
到这一步矩阵的计算都不成问题
d=subs(d,{'r','l','x','y','z','d1','d2','d3'},{31.8,40.2,B(n,1),B(n,2),B(n,3),B(n,4),B(n,5),B(n,6)});
%给参数赋值，其中x y z d1 d2 d3是用B矩阵的第n行对应的列来赋值的
e = norm(d,'fro'); %求范数
f(n,1)=e;%将范数记录
n=n+1;
if n==101;%B数组有100行所以到一百终止
end

最后结果f数列中只有一个数 ，请帮忙
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼2013-04-15 09:14:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

材料廖

木虫 (正式写手)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
加油大楠楠: 金币+10, ★★★★★最佳答案 2013-04-26 10:18:01
引用回帖:
10楼: Originally posted by 加油大楠楠 at 2013-04-23 09:11:04
我昨天晚上研究了一晚晚上  ，都不行说是，我开始是用surf直接画图，把hh所有得数选定了surf（hh(1:34),hh(1:21)),这个命令就能画出来，用pol2cart 我是这么编的for i=1:1:34;
for j=1:1:21;
=pol2cart(i*cos(j* ...

我想你应该是对rad，和theta进行优化吧，对于rad和theta是可以画图的，但是我只能绘制直角坐标的图，柱坐标的我也没想出来，但是直角坐标也应该可以看出来。你在后面添加如下语句：
theta = 0:pi/10:2*pi;
rad = 1:34;
surf(theta,rad,hh);
这样可以得出个以theta为x轴，rad为y轴，hh为z轴的图，你看是不是你要的。
赞一下(1人) 回复此楼
11楼2013-04-23 13:16:50
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

材料廖

木虫 (正式写手)

【答案】应助回帖

感谢参与，应助指数 +1
你这个肯定只有n=1的值呀，都没有循环，你可以将你的这段程序
f(n,1)=e;%将范数记录
n=n+1;
if n==101;%B数组有100行所以到一百终止
end
改成：
for n = 1:100
f(n,1) = e;
end
这样应该就可以了！
赞一下 回复此楼
2楼2013-04-15 12:12:41
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

加油大楠楠

铜虫 (小有名气)
引用回帖:
2楼: Originally posted by 材料廖 at 2013-04-15 12:12:41
你这个肯定只有n=1的值呀，都没有循环，你可以将你的这段程序
f(n,1)=e;%将范数记录
n=n+1;
if n==101;%B数组有100行所以到一百终止
end
改成：
for n = 1:100
f(n,1) = e;
end
这样应该就可以了！

不好意思回复晚了  那你再帮我看看这整段程序   我把上次的程序稍微改了下，就是对于每个l，我要求他xyx为不同值的时候的范数，然后相加，最后把每个l的范数和提取出来，十分感谢了  我运行总是报错
j=1;z=20;r=32;
for l=32:34   %设定范围
i=1;  tezheng=0;     
    for rad=0:l; %这两个for的步长可以根据运算效率适当放大
        for thea=0:pi/10:2*pi;

m=rad*cos(thea);
n=rad*sin(thea);
d11=1/3*(3*r+3*m+sqrt(9*l^2-9*n^2-9*z^2));
d12=1/3*(6*r-3*m+3*n*sqrt(3)+sqrt(36*l^2-27*m^2-18*m*n*sqrt(3)-9*n^2-36*z^2));
d13=1/3*(6*r-3*m-3*n*sqrt(3)+sqrt(36*l^2-27*m^2+18*m*n*sqrt(3)-9*n^2-36*z^2));
B(i,=[m,n,z,d11,d12,d13];
syms x1 y1 z1 d1 d2 d3 l1;
a=(1/l1)*[z1+r-d1 0 0;0 -0.866*y1-0.5*z1+r-d2 0;0 0 0.866*y1-0.5*z1+r-d3];
b=[x1 y1 z1+r-d1;x1 y1-0.866*r+0.866*d2 z1-0.5*r+0.5*d2;x1 y1+0.866*r-0.866*d3 z1-0.5*r+0.5*d3];
c=inv(a);
d=c*b;%d为雅克比矩阵
d=subs(d,{'l1','x1','y1','z1','d1','d2','d3'},{l,B(i,1),B(i,2),B(i,3),B(i,4),B(i,5),B(i,6)});
e = norm(d,'fro'); %2范数
hh(i,1)=1/e;%把这个范数记录在hh中
tezheng=hh(i,1)+tezheng;%相加对于相同l值的对应的范数
   i=i+1;
        end
    end
tezheng1(j,=tezheng/(i-1);%把不同l值对应的范数求和再除以相应的l值  存在数组tezheng1中
j=j+1;
end


错误信息是
clear
??? Error using ==> mupadmex
Error in MuPAD command: Out of memory

Error in ==> sym.horzcat at 25
    y = mupadmex('mllib::horzcat',strs{:});

Error in ==> jacobi at 16
a=(1/l1)*[z1+r-d1 0 0;0 -0.866*y1-0.5*z1+r-d2 0;0 0
0.866*y1-0.5*z1+r-d3];

>>
赞一下 回复此楼
3楼2013-04-18 09:38:21
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

材料廖

木虫 (正式写手)

【答案】应助回帖

我还是不是很明白你要做什么，你这个程序写的我有点看不懂。但有点建议，一：你这个应该没有必要用符号计算，符号计算很慢。二：你这个程序编写的不怎样，基础不是很好，写的有C语言的感觉，可能你的想法你没用程序表达出来，要不你说说你要干啥吧，我看能不能帮你。话说直了点，望不怪。
下面是我没用符号计算程序，速度快了很多，但是还是存在错误，出不了结果，因为我还是没有理解你要干啥。
j=1;z=20;r=32;
for l=32:34   %设定范围
    i=1;  tezheng=0;
    for rad=0:l; %这两个for的步长可以根据运算效率适当放大
        for thea=0:pi/10:2*pi;
            
            m=rad*cos(thea);
            n=rad*sin(thea);
            d11=1/3*(3*r+3*m+sqrt(9*l^2-9*n^2-9*z^2));
            d12=1/3*(6*r-3*m+3*n*sqrt(3)+sqrt(36*l^2-27*m^2-18*m*n*sqrt(3)-9*n^2-36*z^2));
            d13=1/3*(6*r-3*m-3*n*sqrt(3)+sqrt(36*l^2-27*m^2+18*m*n*sqrt(3)-9*n^2-36*z^2));
            B(i,=[m,n,z,d11,d12,d13];
            a=(1/B(i,1))*[B(i,3)+r-B(i,4) 0 0;0 -0.866*B(i,2)-0.5*B(i,3)+r-B(i,5) 0;0 0 0.866*B(i,2)-0.5*B(i,3)+r-B(i,6)];
            b=[B(i,1) B(i,2) B(i,3)+r-B(i,4);B(i,1) B(i,2)-0.866*r+0.866*B(i,5) B(i,3)-0.5*r+0.5*B(i,5);B(i,1) B(i,2)+0.866*r-0.866*B(i,6) B(i,3)-0.5*r+0.5*B(i,6)];
%             c=inv(a);
            d=a\b;%d为雅克比矩阵
            e = norm(d,'fro'); %2范数
            hh(i,1)=1/e;%把这个范数记录在hh中
            tezheng=hh(i,1)+tezheng;%相加对于相同l值的对应的范数
            i=i+1;
        end
    end
    tezheng1(j,=tezheng/(i-1+eps);%把不同l值对应的范数求和再除以相应的l值  存在数组tezheng1中
    j=j+1;
end
赞一下 回复此楼
4楼2013-04-19 18:36:16
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

加油大楠楠

铜虫 (小有名气)
引用回帖:
4楼: Originally posted by 材料廖 at 2013-04-19 18:36:16
我还是不是很明白你要做什么，你这个程序写的我有点看不懂。但有点建议，一：你这个应该没有必要用符号计算，符号计算很慢。二：你这个程序编写的不怎样，基础不是很好，写的有C语言的感觉，可能你的想法你没用程序 ...

看来真是行家  我这是临时学的matlab 才不到一个月  因为我有一个计算要做我跟你说一下  麻烦您帮我看一下，我是做机械的，里面有一个设计参数l，有一个设计指标，是一个矩阵就是上面说的那个雅克比矩阵，我现在要计算对应每个l的雅克比矩阵，但是每个l的雅克比矩阵同时也有变量xyz，d1d2d3的，并不是一个数，他是一个数集，所以你看到我编程用的是i=1;  tezheng=0;     
    for rad=0:l; %这两个for的步长可以根据运算效率适当放大
        for thea=0:pi/10:2*pi;

m=rad*cos(thea);
n=rad*sin(thea);
d11=1/3*(3*r+3*m+sqrt(9*l^2-9*n^2-9*z^2));
d12=1/3*(6*r-3*m+3*n*sqrt(3)+sqrt(36*l^2-27*m^2-18*m*n*sqrt(3)-9*n^2-36*z^2));
d13=1/3*(6*r-3*m-3*n*sqrt(3)+sqrt(36*l^2-27*m^2+18*m*n*sqrt(3)-9*n^2-36*z^2));
说白了我就是想把每个l中不同的xyz，d1d2d3对应的矩阵二范数求出来，就是这样  谢谢了
赞一下 回复此楼
5楼2013-04-21 10:15:17
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

加油大楠楠

铜虫 (小有名气)
引用回帖:
4楼: Originally posted by 材料廖 at 2013-04-19 18:36:16
我还是不是很明白你要做什么，你这个程序写的我有点看不懂。但有点建议，一：你这个应该没有必要用符号计算，符号计算很慢。二：你这个程序编写的不怎样，基础不是很好，写的有C语言的感觉，可能你的想法你没用程序 ...

程序中的吗，m,n就是xy,为了减少计算量我把z固定了是20，其实xyz是空间坐标
赞一下 回复此楼
6楼2013-04-21 10:17:29
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

材料廖

木虫 (正式写手)

【答案】应助回帖

★ ★
jjdg: 金币+2, 辛苦了 2013-04-22 10:54:37
引用回帖:
6楼: Originally posted by 加油大楠楠 at 2013-04-21 10:17:29
程序中的吗，m,n就是xy,为了减少计算量我把z固定了是20，其实xyz是空间坐标...

开始那个由于你在rad中用了0，导致了出现了NAN，我照你的程序修改了下，可以出部分结果，程序如下：
j=1;z=20;r=32;
for l=32:34   %设定范围
    i=1;  tezheng=0;
    for rad=1:l; %这两个for的步长可以根据运算效率适当放大，此处我改成了从1开始
        for thea=0:pi/10:2*pi;
            
            x=rad*cos(thea);  %x
            y=rad*sin(thea);  %y
            d11=1/3*(3*r+3*x+sqrt(9*l^2-9*y^2-9*z^2));
            d12=1/3*(6*r-3*x+3*y*sqrt(3)+sqrt(36*l^2-27*x^2-18*x*y*sqrt(3)-9*y^2-36*z^2));
            d13=1/3*(6*r-3*x-3*y*sqrt(3)+sqrt(36*l^2-27*x^2+18*x*y*sqrt(3)-9*y^2-36*z^2));
            B(i,:)=[x,y,z,d11,d12,d13];
            a=(1/B(i,1))*[B(i,3)+r-B(i,4) 0 0;0 -0.866*B(i,2)-0.5*B(i,3)+...
                r-B(i,5) 0;0 0 0.866*B(i,2)-0.5*B(i,3)+r-B(i,6)];
            b=[B(i,1) B(i,2) B(i,3)+r-B(i,4);B(i,1) B(i,2)-0.866*r+...
                0.866*B(i,5) B(i,3)-0.5*r+0.5*B(i,5);B(i,1) B(i,2)+...
                0.866*r-0.866*B(i,6) B(i,3)-0.5*r+0.5*B(i,6)];
%             c=inv(a);
            d=a\b; %d为雅克比矩阵
            e = norm(d,'fro'); %2范数
            hh(i,1)=1/(eps+e);%把这个范数记录在hh中
            tezheng=hh(i,1)+tezheng;%相加对于相同l值的对应的范数
            i=i+1;
        end
    end
    tezheng1(j,:)=tezheng/(i-1+eps);%把不同l值对应的范数求和再除以相应的l值  存在数组tezheng1中
    j=j+1;
end

根据你说的和我理解的，我自己写了一个，运行速度理论上会快点，也可以把你计算过程中的大部分数据记录下来，程序如下：

%对于l我感觉你要是后面是求和的话，没有必要那样，你的求和在多重循环里面，重复了很多次，我就将l设成了你程序中的最大值，这个你自己再考虑下，看着修改吧。

j=1;r=32;l = 34;  %注意我把l设成了一个常数
for rad = 1:34
    %theta = 0:pi/10:2*pi
    x(rad,:) = rad*cos(0:pi/10:2*pi);  %x
    y(rad,:) = rad*sin(0:pi/10:2*pi);  %y
    z(rad,:) = 20*ones(size(0:pi/10:2*pi));
    d11(rad,:)=1/3*(3*r+3*x(rad,:)+sqrt(9*l^2-9*y(rad,:).^2-9*z(rad,:).^2));
    d12(rad,:)=1/3*(6*r-3*x(rad,:)+3*y(rad,:)*sqrt(3)+sqrt(36*l^2-...
        27*x(rad,:).^2-18*x(rad,:).*y(rad,:)*sqrt(3)-9*y(rad,:).^2-36*z(rad,:).^2));
    d13(rad,:)=1/3*(6*r-3*x(rad,:)-3*y(rad,:)*sqrt(3)+sqrt(36*l^2-...
        27*x(rad,:).^2+18*x(rad,:).*y(rad,:)*sqrt(3)-9*y(rad,:).^2-36*z(rad,:).^2));
    B(rad,1:21,1:6) = [(x(rad,:))',(y(rad,:))',(z(rad,:))',(d11(rad,:))',(d12(rad,:))',(d13(rad,:))'];
    n = length(B(rad,:,1));
    for k = 1:n
        a(:,:,rad,k) = (1/B(rad,k,1)).*[B(rad,k,3)+r-B(rad,k,4) 0 0;0 -0.866.*B(rad,k,2)-0.5.*B(rad,k,3)+...
            r-B(rad,k,5) 0;0 0 0.866.*B(rad,k,2)-0.5.*B(rad,k,3)+r-B(rad,k,6)];
        b(:,:,rad,k) = [B(rad,k,1) B(rad,k,2) B(rad,k,3)+r-B(rad,k,4);B(rad,k,1) B(rad,k,2)-0.866.*r+...
            0.866.*B(rad,k,5) B(rad,k,3)-0.5.*r+0.5.*B(rad,k,5);B(rad,k,1) B(rad,k,2)+...
            0.866.*r-0.866.*B(rad,k,6) B(rad,k,3)-0.5.*r+0.5.*B(rad,k,6)];
        d(:,:,rad,k)=a(:,:,rad,k)\b(:,:,rad,k); %雅克比矩阵
        e(rad,k) = norm(d(:,:,rad,k),'fro');  %二范数
    end
    hh = 1./e;  %同你程序中的hh，只是此处为一个矩阵，即所有你的hh都在这里
end
tezheng = sum((sum(hh)));%所有的hh求和
%到这里，求和都做了，所有的e，hh都有数据在里面，后面我不是很明白，你自己看着修改吧。
赞一下 回复此楼

» 本帖已获得的红花（最新10朵）

加油大楠楠
7楼2013-04-21 18:10:57
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

加油大楠楠

铜虫 (小有名气)
送红花一朵
引用回帖:
7楼: Originally posted by 材料廖 at 2013-04-21 18:10:57
开始那个由于你在rad中用了0，导致了出现了NAN，我照你的程序修改了下，可以出部分结果，程序如下：
j=1;z=20;r=32;
for l=32:34   %设定范围
    i=1;  tezheng=0;
    for rad=1:l; %这两个for的步长可以根据 ...

我大概明白了  再厚脸皮问你一下 这个hh矩阵能不能用极坐标画出来 画成三维的图，就是图中x坐标为极轴长，y坐标为角度，hh中对应的数值为z值
赞一下 回复此楼
8楼2013-04-22 20:32:35
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

材料廖

木虫 (正式写手)
引用回帖:
8楼: Originally posted by 加油大楠楠 at 2013-04-22 20:32:35
我大概明白了  再厚脸皮问你一下 这个hh矩阵能不能用极坐标画出来 画成三维的图，就是图中x坐标为极轴长，y坐标为角度，hh中对应的数值为z值...

在此程序中，x，y均为34*21行的矩阵，而hh是34*21的矩阵，是不能直接绘制三维图的。对于每一重循环，x，y，hh貌似均为向量也是不行的。你自己可以修改看看，要是要用柱坐标绘图，你可以用pol2cart，具体用法可以help pol2cart，将其转换成直角坐标，再在直角坐标中绘制图形。
赞一下 回复此楼
9楼2013-04-22 22:01:14
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

加油大楠楠

铜虫 (小有名气)
引用回帖:
9楼: Originally posted by 材料廖 at 2013-04-22 22:01:14
在此程序中，x，y均为34*21行的矩阵，而hh是34*21的矩阵，是不能直接绘制三维图的。对于每一重循环，x，y，hh貌似均为向量也是不行的。你自己可以修改看看，要是要用柱坐标绘图，你可以用pol2cart，具体用法可以he ...

我昨天晚上研究了一晚晚上  ，都不行说是，我开始是用surf直接画图，把hh所有得数选定了surf（hh(1:34),hh(1:21)),这个命令就能画出来，用pol2cart 我是这么编的for i=1:1:34;
for j=1:1:21;
[x,y]=pol2cart(i*cos(j*pi/10),i*sin(j*pi/10));
end
end
surf(x,y,hh);就不对  能帮我看看吗
赞一下 回复此楼
10楼2013-04-23 09:11:04
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 加油大楠楠 的主题更新
121/2返回列表12下一页
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭