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

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Matlab两个小问题
111/2返回列表12下一页
查看: 1716  |  回复: 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的回帖

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的回帖

sudo

木虫 (正式写手)
还是贴出你的全部程序吧
赞一下 回复此楼
5楼2011-10-13 20:28:33
 已阅   回复此楼   关注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的回帖

zhmdream

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

下面是输入的参数：

输入平台系统参数
太阳光吸收涂层宽度d(m):0.06
玻璃管的直径D(m):0.12
反射镜宽度(m):3
平台系统长度(m):6
太阳光选择吸收涂层的辐射率:0.05
玻璃的辐射率:0.2
反射镜的辐射率:0.2
冷却管背部的辐射率:0.2
冷却管背部的吸收率:0.8
输入周围环境参数
周围环境风速(m/s):1.5
环境空气温度Ta(K)：300
冷却管温度Tc(K):323
赞一下 回复此楼
人生如梦，岁月无情~~~
7楼2011-10-14 10:54:29
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

sudo

木虫 (正式写手)

【答案】应助回帖

★ ★
jjdg(金币+2): 感谢参与 2011-10-15 01:56:01
嗯，这样，其实呢，matlab的solve函数的能力是有限的，虽然一元四次方程一定能解出解析解，但是solve却不一定能算出来

比如
CODE:
solve('x^4+x+c=0','x')

所以对于你的问题中的略显复杂的方程（在matlab看来~），最好的方式是计算数值解

用到的函数是fzero，嗯~
程序如下，你可以自己去体会一下：
CODE:
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('冷却管背部的吸收率:');
d=0.06;
D=0.12;
D_tro=3;
L=6;
E_p=0.05;
E_g=0.2;
E_tro=0.2;
E_rear=0.2;
a_rear=0.8;

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

%常数
%斯蒂芬－波尔兹曼常数
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);

%求数值解
tp = 450:650; %因为上面f里面的Tp还有用，所以这里不能同名
tg = zeros(1, length(tp));

for i = 1:length(tp)
    tf = subs(f, Tp, tp(i)); %给Tp代入值
    fHandle = inline(char(tf), 'Tg'); %将tf视为函数句柄，'Tg'为自变量
   
    %在Tg=400附近的实数域找解，如果找不到会返回NaN
    tg(i) = fzero(fHandle, 400);
   
    if tg(i)<=400
        tg(i) = 0; %如果没找到大于400的解，那么就置为0表示没找到
    end
end

plot(tp, tg);

之后画出的图是这个样子：
赞一下(1人) 回复此楼
8楼2011-10-14 13:02:07
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

sudo

木虫 (正式写手)
★
jjdg(金币+1): 好建议 2011-10-15 01:56:10
要善用help命令和google啊...

上面这个程序里面，需要注意的手法只有一个：

如何把一个符号表达式，转化为函数句柄（因为MATLAB里面很多解方程/优化函数里面都有一个“函数句柄”的参数）
赞一下(1人) 回复此楼
9楼2011-10-14 13:07:46
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

我有风衣

银虫 (正式写手)
真是不错！编程新人来学习一下！谢谢！
赞一下 回复此楼
南理工--FPGA应用
10楼2011-12-01 09:50:56
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 zhmdream 的主题更新
111/2返回列表12下一页
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭