24小时热门版块排行榜    

查看: 2031  |  回复: 6

watertxf

铁虫 (初入文坛)

[求助] matlab解四次方程的问题 已有1人参与

matlab解四次方程的问题,我写的程序如下:要求在w=0.9:0.01:1.1的范围内求x的最大值*R随gamma=1.1:0.5:10的变化,大家能帮忙看一下出什么问题了吗?
非常感谢!
clear
clc
syms x
w=0.9:0.01:1.1;%21个
gamma=1.1:0.5:10;%18个
c=3*10^8;
wpb=0.2;
R=2.3*10^-2;


for j=1:18 %gamma
    v(j)=c*sqrt(1-1/gamma(j)^2);
    for i=1:21 %w
    y(i)=1+0.2^2/(1-w(i)^2+2i*0.015*w(i));
    %wpb(j)=sqrt((1.602*10^-19)^2*10^9/(9.109*10^-31*8.85*10^-12*gamma(j)^3));
    f(i,j)=(y(i)*(w(i)*2*pi*24*10^9)^2/c^2-x^2)*(1-(wpb*2*pi*24*10^9)^2/(((w(i)*2*pi*24*10^9)-x*v(j))^2*y(i)))-2.4048^2/R^2;
    k=solve(f);  %解出来k有四个根,
    k2=imag(k)*R;%需要k2的第四个根的值
    k2=double(k2);
    %m1(i)=k1(1,;
    %m2(i)=k1(2,;
    %m3(i)=k1(3,;
    %m4(i)=k1(4,;
    %n1(i)=k2(1,;
    %n2(i)=k2(2,;
    %n3(i)=k2(3,;
    %n4(i)=k2(4,;
    end
   k3(j)=max(k2);
end
%d=max(d)

semilogy(gamma,k3,'b --')
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hppdyx

木虫 (知名作家)

【答案】应助回帖

我没有试,也不知道出的什么错误,不过solve求解都是把变量当作sys来看,所以求解出来的解是表达式而不是数值,因此可以试着加上eval函数,看看行不行

» 本帖已获得的红花(最新10朵)

不以风骚惊天下,但求淫荡动世人
2楼2013-12-06 13:09:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

watertxf

铁虫 (初入文坛)

送红花一朵
引用回帖:
2楼: Originally posted by hppdyx at 2013-12-06 13:09:21
我没有试,也不知道出的什么错误,不过solve求解都是把变量当作sys来看,所以求解出来的解是表达式而不是数值,因此可以试着加上eval函数,看看行不行

非常感谢您的帮助!我对matlab不是太熟悉,能不能麻烦您告诉我怎样加上eval函数求解?不胜感激!
3楼2013-12-10 14:32:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hppdyx

木虫 (知名作家)

引用回帖:
3楼: Originally posted by watertxf at 2013-12-10 14:32:57
非常感谢您的帮助!我对matlab不是太熟悉,能不能麻烦您告诉我怎样加上eval函数求解?不胜感激!...

首先,你这个程序里面有很多小细节上的错误,还有一些错误需要你解释一下才行。我把你的程序改了一下,调试之后可以了,不过不知道符不符合你的要求。
CODE:
function k = res2

syms x
w = 0.9 : 0.01 : 1.1; %21个
gamma = 1.1 : 0.5 : 10;   %18个
c = 3 * 10^8;

R = 2.3 * 10^-2;

n1 = length(w);
n2 = length(gamma);
v = zeros(n2);
y = zeros(n1);
wpb = zeros(n2);
f = sym(zeros(n1, n2));
k = zeros(n1 * n2, 1);

for j = 1 : 18        %gamma
    v(j) = c * sqrt(1 - 1 ./ gamma(j).^2);
    for i = 1 : 21        %w
    y(i) = 1 + 0.2^2 ./ (1 - w(i).^2 + 2 * i * 0.015 * w(i));
    wpb(j) = sqrt((1.602 * 10^-19)^2 * 10^9 ./ (9.109 * 10^-31 * 8.85 * 10^-12 * gamma(j)^3));
    f(i, j) = (y(i) * (w(i) * 2 * pi * 24 * 10^9)^2 / c^2 - x^2) * (1 - (wpb(j) * 2 * pi * 24 * 10^9)^2 / (((w(i) * 2 * pi * 24 * 10^9) - x * v(j))^2 * y(i))) - 2.4048^2 / R^2;
    k((i-1) * n2 + j) = max(sqrt(eval(solve(f(i, j)))));  
    end
end
semilogy(gamma, k(1:18), 'b--');
figure;
plot(gamma, k(1:18), 'r');

结果如图所示:
matlab解四次方程的问题
matlab解四次方程的问题-1
不以风骚惊天下,但求淫荡动世人
4楼2013-12-12 20:26:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hppdyx

木虫 (知名作家)

引用回帖:
4楼: Originally posted by hppdyx at 2013-12-12 20:26:07
首先,你这个程序里面有很多小细节上的错误,还有一些错误需要你解释一下才行。我把你的程序改了一下,调试之后可以了,不过不知道符不符合你的要求。
function k = res2

syms x
w = 0.9 : 0.01 : 1.1; %21个 ...

不好意思,倒数第六行的sqrt改为abs。。。。图形的走势是一致的
不以风骚惊天下,但求淫荡动世人
5楼2013-12-12 20:43:43
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hppdyx

木虫 (知名作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
watertxf: 金币+20, ★★★★★最佳答案, 感谢! 2013-12-13 15:21:57
引用回帖:
4楼: Originally posted by hppdyx at 2013-12-12 20:26:07
首先,你这个程序里面有很多小细节上的错误,还有一些错误需要你解释一下才行。我把你的程序改了一下,调试之后可以了,不过不知道符不符合你的要求。
function k = res2

syms x
w = 0.9 : 0.01 : 1.1; %21个 ...

把sqrt改为abs,然后又把点加密了一些,程序和结果如图(一个纵坐标是对数坐标,另一个纵坐标不是对数坐标):
CODE:
function k = res2

syms x
w = 0.9 : 0.01 : 1.1;
gamma = 1.1 : 0.2 : 10;  
c = 3 * 10^8;

R = 2.3 * 10^-2;

n1 = length(w);
n2 = length(gamma);
v = zeros(n2);
y = zeros(n1);
wpb = zeros(n2);
f = sym(zeros(n1, n2));
k = zeros(n1 * n2, 1);

for j = 1 : n2       
    v(j) = c * sqrt(1 - 1 ./ gamma(j).^2);
    for i = 1 : n1       
    y(i) = 1 + 0.2^2 ./ (1 - w(i).^2 + 2 * i * 0.015 * w(i));
    wpb(j) = sqrt((1.602 * 10^-19)^2 * 10^9 ./ (9.109 * 10^-31 * 8.85 * 10^-12 * gamma(j)^3));
    f(i, j) = (y(i) * (w(i) * 2 * pi * 24 * 10^9)^2 / c^2 - x^2) * (1 - (wpb(j) * 2 * pi * 24 * 10^9)^2 / (((w(i) * 2 * pi * 24 * 10^9) - x * v(j))^2 * y(i))) - 2.4048^2 / R^2;
    k((i-1) * n2 + j) = max(abs(eval(solve(f(i, j)))));  
    end
end
semilogy(gamma, k(1:n2), 'b--');
figure;
plot(gamma, k(1:n2), 'r');

matlab解四次方程的问题-2
matlab解四次方程的问题-3
不以风骚惊天下,但求淫荡动世人
6楼2013-12-12 20:52:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

watertxf

铁虫 (初入文坛)

引用回帖:
6楼: Originally posted by hppdyx at 2013-12-12 20:52:09
把sqrt改为abs,然后又把点加密了一些,程序和结果如图(一个纵坐标是对数坐标,另一个纵坐标不是对数坐标):
function k = res2

syms x
w = 0.9 : 0.01 : 1.1;
gamma = 1.1 : 0.2 : 10;  
c = 3 * 10^8 ...

太感谢了!这个趋势是正确的,只是数量级不对,我在看一下。非常感谢您的帮助!
7楼2013-12-13 15:22:35
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 watertxf 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 面上函评意见出来了,像什么等级? 20+4 Tsingking1 2026-08-27 17/850 2026-09-01 19:51 by 超级无敌华子
[文学芳草园] 梦想 +5 myrtle 2026-08-26 7/350 2026-09-01 15:18 by myrtle
[基金申请] 学科评审组评审是指会评吗? +4 瞬息宇宙 2026-08-31 4/200 2026-09-01 14:58 by jiaoxg
[基金申请] 为什么资助数各大高校都创新高,自己申请怎么就这么难 +13 Kittylucky 2026-08-27 14/700 2026-09-01 11:06 by feng6531
[论文投稿] 小白求助 投论文要求的highlights应该如何写 5+3 l1963982152 2026-08-29 4/200 2026-09-01 09:04 by 北京莱茵编辑
[基金申请] 国社科又开始会评了,不知道这次命运如何 +7 雨打竹帘 2026-08-30 11/550 2026-08-31 23:16 by hittle2008
[基金申请] 投票:  有多少人是今天查系统知道结果的? +17 爱看书的可乐 2026-08-26 19/950 2026-08-31 21:30 by xiangy672
[基金申请] 面上意见出来了 +12 黄鸟于飞Chao 2026-08-29 23/1150 2026-08-31 18:57 by 黄鸟于飞Chao
[基金申请] 中青基了要发朋友圈吗? +7 349506619 2026-08-28 7/350 2026-08-31 13:39 by 冼亮淀粉酶
[基金申请] 基金不中,共勉 +12 eulota 2026-08-26 12/600 2026-08-31 08:42 by ZJTJZ
[基金申请] 有没有仍没收到信息的 +7 德尚中行 2026-08-27 8/400 2026-08-30 20:52 by purplejack
[基金申请] 我就是申请一个面上项目而已,这评审意见是按照杰青的条件评的吧? +6 gouxfjh 2026-08-28 11/550 2026-08-30 07:57 by gouxfjh
[基金申请] 国自然面上复盘~欢迎讨论 (金币+15) +15 晴天加油 2026-08-26 16/800 2026-08-29 18:28 by symmetry
[基金申请] 国自然评审意见 +13 wangmingqi 2026-08-28 19/950 2026-08-29 10:22 by Poppy1104
[基金申请] 怎么查啊 +6 huang1991js 2026-08-26 6/300 2026-08-28 08:42 by winsaint
[基金申请] 看板上这么多中的,有点像50人群里49个人都是骗子的那种感觉…… +5 a089 2026-08-26 6/300 2026-08-27 14:05 by jonewore
[基金申请] 国际合作可查了,中了面上 (EPI+1)(金币+50) +18 Ldrop2023 2026-08-26 18/900 2026-08-26 11:15 by cmrandy
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
[基金申请] 国合可查了 +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
信息提示
请填处理意见