24小时热门版块排行榜    

查看: 1385  |  回复: 14
当前主题已经存档。
本帖产生 1 个 仿真EPI ,点击这里进行查看

feibao

银虫 (小有名气)

真的很谢谢你!

如开篇程序所言在MATLAB10版本中对行列式总体求导所得的解析式也是能够较快求出,但在赋值时迟迟没有进展;但分开求导,哪怕所得解析式仍汇合在一起 ,最后对总解析式赋值,却仍能较快求出数值。不知何故???

两种方法对最后赋值的影响难道机理上有什么玄机?

clear all
format long
syms u  v

m=1;
s=1+0.1;
a=0.98;
b=0.99;
n1=1.48;
n2=-a*n1;
n3=b*n1;
r=(n1^2-n2^2)/(n1^2-n3^2);

w1=v^2*r-u^2;
w=v^2-u^2;
ne=sqrt((1-u^2/v^2)*(n1^2-n3^2)+n3^2);
F1=-besselj(m+1,u)+m./u.*besselj(m,u);
F=F1./(u.*besselj(m,u));
M1=besseli(m+1,w1)+m./w1.*besseli(m,w1);
M=M1./(w1.*besseli(m,w1));
N1=-besselk(m+1,w1)+m./w1.*besselk(m,w1);

N=N1./(w1.*besselk(m,w1));
Q=m*(1./u.^2+1./w1.^2);
R=m*(1./w1.^2-1./w.^2)/s;
S1=besseli(m+1,w1*s)+m./w1./s.*besseli(m,w1*s);

S=S1./(w1.*besseli(m,w1*s));
T1=-besselk(m+1,w*s)+m./w./s.*besselk(m,w*s);
T=T1./(w.*besselk(m,w*s));
X1=-besselk(m+1,w1*s)+m./w1./s.*besselk(m,w1*s);

X=X1./(w1.*besselk(m,w1*s));
I1m=besseli(m,w1);
K1m=besselk(m,w1);
I2m=besseli(m,w1*s);
K2m=besselk(m,w1*s);
x1=Q.*ne.*I1m;
x2=Q.*ne.*K1m;
x3=I1m.*(F-M);
x4=K1m.*(F-N);
x5=R*ne.*I2m;
x6=R.*ne.*K2m;
x7=I2m.*(-S-T);
x8=K2m.*(-X-T);
x9=I1m.*(n1^2.*F-n2^2.*M);
x10=K1m.*(n1^2.*F-n2^2.*N);
x11=Q.*ne.*I1m;
x12=Q.*ne.*K1m;
x13=I2m.*(-S.*n2^2-n3^2.*T);
x14=K2m.*(-n2^2.*X-n3^2.*T);
x15=R.*ne.*I2m;
x16=R.*ne.*K2m;
% f=det([x1,x2,x3,x4;
%        x5,x6,x7,x8;
%        x9,x10,x11,x12;
%        x13,x14,x15,x16]);

%%PSL@CSU
%%QQ:547423688
%%Email:anyuezhiji@qq.com
%%Edit @ 2010.4.17

%拆开算  
detf1=x1*x6*x11*x16;
detf2=-x1*x6*x12*x15;
detf3=-x1*x10*x7*x16;
detf4=x1*x10*x8*x15;
detf5=x1*x14*x7*x12;
detf6=-x1*x14*x8*x11;
detf7=-x5*x2*x11*x16;
detf8=x5*x2*x12*x15;
detf9=x5*x10*x3*x16;
detf10=-x5*x10*x4*x15;
detf11=-x5*x14*x3*x12;
detf12=x5*x14*x4*x11;
detf13=x9*x2*x7*x16;
detf14=-x9*x2*x8*x15;
detf15=-x9*x6*x3*x16;
detf16=x9*x6*x4*x15;
detf17=x9*x14*x3*x8;
detf18=-x9*x14*x4*x7;
detf19=-x13*x2*x7*x12;
detf20=x13*x2*x8*x11;
detf21=x13*x6*x3*x12;
detf22=-x13*x6*x4*x11;
detf23=-x13*x10*x3*x8;
detf24=x13*x10*x4*x7;

for i=1:24
eval(['fv',num2str(i),'=diff(detf',num2str(i),',v);']);
end
for i=1:24
eval(['fu',num2str(i),'=diff(detf',num2str(i),',u);']);
end

%这里开始读入数据
% v=6.4414247e-001   
%u=6.4091366e-001  
%q=load('e:\0.1.txt');

%for i=1:length(q(:,1))
% u=q(i,3)
   % v=q(i,1)

%下面开始求值
ffv=0;
for i=1:24
    eval(['ffv=ffv+eval(fv',num2str(i),');'])
end

ffu=0;
for i=1:24
    eval(['ffu=ffu+eval(fu',num2str(i),');'])
end

ff1=-ffv/ffu;
ff2=1+(u/v)^2*(1-2*v/u*ff1);
q=load('e:\0.1.txt');
for ii=1:length(q(:,1))

v=q(ii,1)
u=q(ii,3)
ff=eval(ff2)

hold on

plot(v,ff)
end
11楼2010-04-17 14:28:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

feibao

银虫 (小有名气)

引用回帖:
Originally posted by anyuezhiji at 2010-04-17 13:32:24:
用下面的代码测试了下



分开算再相加直接求导得出的结果是一样的

真的很谢谢你!

如开篇程序所言在MATLAB10版本中对行列式总体求导所得的解析式也是能够较快求出,但在赋值时迟迟没有进展;但分开求导(如您所言),哪怕所得解析式仍汇合在一起 ,最后对总解析式赋值,却仍能较快求出数值。不知何故???

两种方法对最后赋值的影响难道机理上有什么玄机?
12楼2010-04-17 14:34:05
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

anyuezhiji

银虫 (正式写手)

星空行者

feibao(金币+3): 2010-04-18 12:16
feibao(金币+2): 2010-04-21 13:58
引用回帖:
Originally posted by feibao at 2010-04-17 14:28:07:
真的很谢谢你!

如开篇程序所言在MATLAB10版本中对行列式总体求导所得的解析式也是能够较快求出,但在赋值时迟迟没有进展;但分开求导,哪怕所得解析式仍汇合在一起 ,最后对总解析式赋值,却仍能较快求出数值 ...

f以及求导后的fv,fu在内存中都是一个变量的形式存储了一串很长的表达式
要求是存在连续的内存中 也许牵扯到在跨段和跨页的问题 在内存装载变量和赋值计算都要消耗内存
分开后可以分成几部分存 不要求是连续内存
计算也是分开的 这样算完一部分就能清空一部分 然后再加载后面的计算
实际对最大连续内存的需求不到原来1/24

单个变量所占内存太大处理起来可能出错


上面我做的处理没有把解析式加到一起再赋值
而是把每个解析式赋值
再把赋值的结果相加

[ Last edited by anyuezhiji on 2010-4-17 at 14:45 ]
暗月下没有留下风的痕迹,但它已经寂然飘逝。。By<暗月之寂>:tiger38:
13楼2010-04-17 14:37:05
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

anyuezhiji

银虫 (正式写手)

星空行者

★ ★
nono2009(金币+2, 仿真EPI+1):综合本帖交流。 2010-04-18 11:17
feibao(金币+5): 2010-04-18 12:20
赋值那里可以改为,这样对内存的要求可能少些:
引用回帖:
%这里开始读入数据
q=load('e:\0.1.txt');
for ii=1:length(q(:,1))

v=q(ii,1)
u=q(ii,3)

%下面开始求值
ffv=0;
for i=1:24
    eval(['ffv=ffv+eval(fv',num2str(i),');'])
end

ffu=0;
for i=1:24
    eval(['ffu=ffu+eval(fu',num2str(i),');'])
end

ff1=-ffv/ffu;
ff2=1+(u/v)^2*(1-2*v/u*ff1);

ff=eval(ff2)

hold on

plot(v,ff)
end

引用回帖:
%这里开始读入数据
q=load('e:\0.1.txt');

ff=zeros(length(q(:,1)));
for ii=1:length(q(:,1))

v=q(ii,1)
u=q(ii,3)

%下面开始求值
ffv=0;
for i=1:24
    eval(['ffv=ffv+eval(fv',num2str(i),');'])
end

ffu=0;
for i=1:24
    eval(['ffu=ffu+eval(fu',num2str(i),');'])
end

ff1=-ffv/ffu;
ff2=1+(u/v)^2*(1-2*v/u*ff1);

ff(ii)=eval(ff2);

end

v=q(:,1);

plot(v,ff);


[ Last edited by anyuezhiji on 2010-4-17 at 15:03 ]
暗月下没有留下风的痕迹,但它已经寂然飘逝。。By<暗月之寂>:tiger38:
14楼2010-04-17 14:56:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

liuyalei820

金虫 (正式写手)


anyuezhiji(金币+1):感谢回帖交流 2010-04-21 14:57
分开算再相加直接求导得出的结果是一样的
15楼2010-04-21 11:07:02
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 feibao 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 面上意见出来了 +6 黄鸟于飞Chao 2026-08-29 10/500 2026-08-29 15:50 by 微不足道士
[教师之家] 导师吐槽:我怎么摊上了这么个极品研究生! +9 苏东坡二世 2026-08-23 9/450 2026-08-29 14:41 by hustersqt
[博后之家] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +4 ASdOkHsho7FD 2026-08-28 5/250 2026-08-29 14:24 by jCd0dEvKHShX
[论文投稿] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +3 ASdOkHsho7FD 2026-08-28 4/200 2026-08-29 14:13 by jCd0dEvKHShX
[硕博家园] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +5 ASdOkHsho7FD 2026-08-28 9/450 2026-08-29 14:03 by jCd0dEvKHShX
[论文投稿] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +3 ASdOkHsho7FD 2026-08-28 4/200 2026-08-29 13:51 by jCd0dEvKHShX
[硕博家园] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +3 ASdOkHsho7FD 2026-08-28 9/450 2026-08-29 13:41 by jCd0dEvKHShX
[考研] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +6 ASdOkHsho7FD 2026-08-28 9/450 2026-08-29 11:18 by jCd0dEvKHShX
[基金申请] 国自然面上复盘~欢迎讨论 (金币+15) +14 晴天加油 2026-08-26 15/750 2026-08-29 11:06 by biily
[基金申请] 国自然评审意见 +13 wangmingqi 2026-08-28 19/950 2026-08-29 10:22 by Poppy1104
[教师之家] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +4 ASdOkHsho7FD 2026-08-28 7/350 2026-08-29 09:26 by G6APbkg8SA6w
[基金申请] 有没有仍没收到信息的 +3 德尚中行 2026-08-27 3/150 2026-08-28 15:24 by rongshuxia
[基金申请] 申请删除本帖 +6 lyz123lyz 2026-08-27 7/350 2026-08-27 17:31 by 宁静致远sy
[文学芳草园] 梦想 +3 myrtle 2026-08-26 3/150 2026-08-27 10:01 by angelyueyi
[基金申请] 为什么 国际(地区)合作与交流项目 没有放榜? 10+3 majunge000 2026-08-26 11/550 2026-08-27 08:42 by 北京莱茵编辑
[基金申请] 2026年的国家社科基金项目通讯评审的新规则与新动向、新挑战 +7 process2012 2026-08-23 10/500 2026-08-26 19:23 by hmhminy
[基金申请] 出来了 +9 trojank 2026-08-26 9/450 2026-08-26 14:25 by 宝贝虫子
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 国合可查了 +3 paperzjh 2026-08-26 3/150 2026-08-26 10:41 by LemmonTr
[基金申请] 牛来!米来!面来! +8 beefly 2026-08-26 8/400 2026-08-26 08:37 by xuzhipiao
信息提示
请填处理意见