请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 3613|回复: 3

[讨论] 求~帮忙编一个一元三次方程计算的matlab程序

[复制链接]
smile123 该用户已被删除
发表于 2019-1-7 21:02:58 | 显示全部楼层 |阅读模式
提示: 作者被禁止或删除 内容自动屏蔽
密码修改失败请联系微信:mofangbao
smile123 该用户已被删除
 楼主| 发表于 2019-1-7 21:03:53 | 显示全部楼层
提示: 作者被禁止或删除 内容自动屏蔽
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-1-7 21:35:54 | 显示全部楼层
function x = Solve3Polynomial(a, b, c, d)
% 范盛金. 一元三次方程的新求根公式与新判别法[J]. 海南师范学院学报,1989,2(2):91-98.
A = b*b - 3*a*c;   if abs(A) < 1e-14;    A = 0;  end
B = b*c - 9*a*d;   if abs(B) < 1e-14;    B = 0;  end
C = c*c - 3*b*d;   if abs(C) < 1e-14;    C = 0;  end  
DET = B*B - 4*A*C; if abs(DET) < 1e-14;  DET = 0;  end  
if (A == 0) && (B == 0)
    x1 = -c/b;      x2 = x1 ;    x3 = x1;
end
if DET > 0
    Y1 = A*b + 1.5*a*(-B + sqrt(DET));
    Y2 = A*b + 1.5*a*(-B - sqrt(DET));
    y1 = nthroot(Y1,3);  y2 = nthroot(Y2,3);
    x1 = (-b-y1-y2)/(3*a);
    vec1 = (-b + 0.5*(y1 + y2))/(3*a);  
    vec2 = 0.5*sqrt(3)*(y1 - y2)/(3*a);
    x2 = complex(vec1, vec2);
    x3 = complex(vec1, -vec2);
    clear Y1 Y2 y1 y2 vec1 vec2;
end
if DET == 0 && (A ~= 0) && (B ~= 0)
    K = (b*c-9*a*d)/(b*b - 3*a*c); K = round(K,14);
    x1 = -b/a + K;   x2 = -0.5*K;   x3 = x2;
end
if DET < 0
    sqA = sqrt(A);
    T = (A*b - 1.5*a*B)/(A*sqA);
    theta = acos(T);
    csth  = cos(theta/3);
    sn3th = sqrt(3)*sin(theta/3);
    x1 = (-b - 2*sqA*csth)/(3*a);
    x2 = (-b + sqA*(csth + sn3th))/(3*a);
    x3 = (-b + sqA*(csth - sn3th))/(3*a);
    clear sqA T theta csth sn3th;
end
x = [x1;  x2;  x3];
end
密码修改失败请联系微信:mofangbao
smile123 该用户已被删除
 楼主| 发表于 2019-1-8 10:42:13 | 显示全部楼层
提示: 作者被禁止或删除 内容自动屏蔽
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表