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

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Matlab两个小问题
南方科技大学公共卫生及应急管理学院2026级博士研究生招生报考通知（长期有效）
51/1返回列表
查看: 1721  |  回复: 10
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
本帖产生 1 个 程序强帖 ，点击这里进行查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

zhmdream

木虫 (正式写手)

[求助] Matlab两个小问题

您好，请教两个问题：
（1）我用Matlab解一个一元4次方程，得到4个解。怎么才能让其输出我想要的实数解，如第一个解（或者解的条件满足：y>400)。

y =
  1.0e+002 *
   4.3541         
   2.9607         
  -3.6574 + 5.2191i
  -3.6574 - 5.2191i

(2)第2个问题是在第1个问题基础上的。我的目的是为了画出x与y的关系图，所以我必须得对连续的x进行求解，然后求出y，作图。以一个简单例子吧：

x=1:1:10; %x取1到10
syms y;
a=y+1;
b=y-1;

f=a*b+a+b*x+y+6*x+2

%解方程f=0
equation = [char(f) ' = ' num2str(0,9)];
A=solve(equation, 'y')

double(A)

plot(x,y)

但是运行出错。有没有其他方法可以解决？谢谢！
回复此楼

» 收录本帖的淘帖专辑推荐

学海无涯！ matlab典型案例及小技巧 大开眼界

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
人生如梦，岁月无情~~~
1楼2011-10-11 20:24:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhmdream

木虫 (正式写手)
引用回帖:
5楼: Originally posted by sudo at 2011-10-13 20:28:33:
还是贴出你的全部程序吧

你好，以下是我的程序：

clear
clc

%输入平台系统参数
display('输入平台系统参数');
d=input('太阳光吸收涂层宽度d(m):');
D=input('玻璃管的直径D(m):');
D_tro=input('反射镜宽度(m):');
L=input('平台系统长度(m):');
E_p=input('太阳光选择吸收涂层的辐射率:');
E_g=input('玻璃的辐射率:');
E_tro=input('反射镜的辐射率:');
E_rear=input('冷却管背部的辐射率:');
a_rear=input('冷却管背部的吸收率:');


%输入周围环境参数
display('输入周围环境参数')
v=input('周围环境风速(m/s):');
Ta=input('环境空气温度Ta(K)：');
Tc=input('冷却管温度Tc(K):');


%常数
%斯蒂芬－波尔兹曼常数
sigma=5.667e-8;  % W/(m2.K4)



%定义玻璃管温度Tg为变量
syms Tg Tp


%%%太阳能热电系统热损耗系数如下%%%
%（1）太阳光选择吸收涂层与玻璃管下部之间辐射：
hr_sg=sigma*(Tp^2+Tg^2)*(Tp+Tg)/(1/E_p+1/E_g-1);
Ap=d*L;
R1=1/(hr_sg*Ap);
%（2）玻璃管上部与冷却管之间的辐射：
hr_rg=sigma*(Tc^2+Tg^2)*(Tc+Tg)/(E_rear+E_g-1);
R2=1/(hr_rg*Ap);
%（3）玻璃管下部与反射镜之间辐射：
hr_gtro=sigma*(Tg^2+Ta^2)*(Tg+Ta)/(E_g+E_tro-1);
Ag=pi*D*L;
R3=1/(hr_gtro*0.5*Ag);
%（4）玻璃管下部与反射镜之间空气对流：
hw=5.7+3.8*v  % W/(m2.K)
R4=1/(hw*0.5*Ag);
%(5)反射镜与周围环境空气对流：
%略
%（6）玻璃管上部对天空的辐射：
hr_gsky=sigma*(Tg^2+Ta^2)*(Tg+Ta)*E_g;
R6=1/(hr_gsky*0.5*Ag);
%（7）反射镜对天空的辐射：
%略
%（8）玻璃管上部与周围环境空气对流：
R8=1/(hw*0.5*Ag);

%求解Tg
f=(Tp-Tg)/R1+(Tc-Tg)/R2-(Tg-Ta)*(1/R3+1/R4+1/R6+1/R8);
equation = [char(f) ' =0 '];
A=solve(equation, 'Tg')

Tp=450:1:650;

for i=1:length(Tp)
    Tg1(i) = subs(A(1), Tp(i));
    Tg2(i) = subs(A(2), Tp(i));
    Tg3(i) = subs(A(3), Tp(i));
    Tg4(i) = subs(A(4), Tp(i));
end


for j=1:length(Tg)
    i=1:4
    if isreal(Tg(i)(j)) && Tg(i)(j)>400 %判断y(i)是否是实数，并且是否大于400
        fprintf(1, '%.9f\n', Tg(i)(j));
    end
end

plot(Tp, Tg)
赞一下 回复此楼
人生如梦，岁月无情~~~
6楼2011-10-14 10:50:55
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 11 个回答

sudo

木虫 (正式写手)

【答案】应助回帖

★ ★ ★ ★
xzhdty(金币+2): 欢迎常来程序语言看看 2011-10-12 15:31:47
微尘、梦想(金币+2, 程序强帖+1): sudo总是很热心的，祝福哦！ 2011-10-12 20:41:46
zhmdream(金币+40): 谢谢啊，我试试 2011-10-12 20:45:04
zhmdream(金币+10): 谢谢了，呵呵 2011-10-14 22:14:37
你不用每次都放那么多分上来的...

(1)使用isreal函数和循环
CODE:
for i=1:length(y)
    if isreal(y(i)) && y(i)>400 %判断y(i)是否是实数，并且是否大于400
        fprintf(1, '%.9f\n', y(i));
    end
end

(2)
CODE:
syms x y; %注意用到的符号都要声明一下
a=y+1;
b=y-1;

f=a*b+a+b*x+y+6*x+2;

%解方程f=0
equation = [char(f) ' = 0'];
A=solve(equation, 'y');

%上面解得的结果是
%A =
%  - x/2 - (x^2 - 16*x - 4)^(1/2)/2 - 1
%   (x^2 - 16*x - 4)^(1/2)/2 - x/2 - 1

%对于符号式子，使用subs函数代入数值

%注意上面的方程因为有开方，所以不一定能解出实数解，这里直接取一段有实数解的
x = 20:80;
y1 = zeros(1, length(x)); %习惯上循环前先给定存储空间
y2 = zeros(1, length(x));
for i=1:length(x)
    y1(i) = subs(A(1), x(i));
    y2(i) = subs(A(2), x(i));
end

plot(x, y1, 'b', x, y2, 'r');
赞一下(2人) 回复此楼
2楼2011-10-12 10:12:19
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

sudo

木虫 (正式写手)
★ ★
微尘、梦想(金币+2): 谢谢参与应助~ 2011-10-12 20:42:05
其实第二个问题里面，算出A之后，直接
CODE:
x = 20:80;
Y= subs(A, x);
plot(x, Y(1, :), 'b', x, Y(2, :), 'r');

就可以的...

不过这么写的话觉得你可能不容易理解
赞一下(1人) 回复此楼
3楼2011-10-12 10:16:33
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhmdream

木虫 (正式写手)
引用回帖:
2楼: Originally posted by sudo at 2011-10-12 10:12:19:
你不用每次都放那么多分上来的...

(1)使用isreal函数和循环
[code]
for i=1:length(y)
    if isreal(y(i)) && y(i)>400 %判断y(i)是否是实数，并且是否大于400
        fprintf(1, '%.9f ...

呵呵，谢谢帮忙。问题解决差不多了，但是还有些小问题：

问题（2）的方程当时是为了方便，随便设的。我要解的是一元四次方程，解的结果有两个实解和两个虚解。舍去虚解，并只有实解>400，才是我所要求的解，并用于作图。下面我的设计与思路：

%%%之前程序省略
%求解Tg
A=solve(equation, 'Tg');   
Tp=450:1:453;     %变量Tp取值，并将Tp代入，求出Tg
for i=1:length(Tp)
    Tg1(i) = subs(A(1), Tp(i));
    Tg2(i) = subs(A(2), Tp(i));
    Tg3(i) = subs(A(3), Tp(i));
    Tg4(i) = subs(A(4), Tp(i));
end
%得到两实解与两虚解

%现在对这两实解与两虚解，进行比较。大于400的实数解才是我最终想要的
for i=1:4
     j=1:length(Tg)
    if isreal(Tg(i)(j)) && Tg(i)(j)>400 %判断Tg(i)是否是实数，并且是否大于400
        fprintf(1, '%.9f\n', Tg(i)(j));
    end
end

plot(Tp, Tg) %对Tp与Tg作图

我运行时，到了第2部分就出错。谢谢！
赞一下 回复此楼
人生如梦，岁月无情~~~
4楼2011-10-13 20:23:16
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 11 个回答
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭