ค้นหาคำตอบทั้งหมดของระบบสมการไม่เชิงเส้นด้วย MuPAD

คำถามของฉันคือมีวิธีที่ดีในการใช้ฟังก์ชัน MuPAD ในสคริปต์ Matlab หรือไม่ เบื้องหลังคือฉันมีปัญหาที่ต้องค้นหาคำตอบทั้งหมดของชุดสมการไม่เชิงเส้น วิธีแก้ปัญหาก่อนหน้านี้คือการใช้ solve ใน Matlab ซึ่งใช้ได้กับการจำลองบางส่วนของฉัน (เช่น ชุดอินพุต T บางชุด) แต่ก็ไม่เสมอไป ดังนั้นฉันจึงใช้ MuPAD แทนด้วยวิธีต่อไปนี้:

function ut1 = testMupadSolver(T)
% # Input T should be a vector of 15 elements

mupadCommand = ['numeric::polysysroots({' eq1(T) ' = 0,' ...
    eq2(T) '= 0},[u, v])'];
allSolutions = evalin(symengine, mupadCommand);
ut1 = allSolutions;

end

function strEq = eq1(T)
sT = @(x) ['(' num2str(T(x)) ')'];

strEq = [ '-' sT(13) '*u^4 + (4*' sT(15) '-2*' sT(10) '-' sT(11) '*v)*u^3 + (3*' ...
    sT(13) '-3*' sT(6) '+v*(3*' sT(14) '-2*' sT(7) ')-' sT(8) '*v^2)*u^2 + (2*' ...
    sT(10) '-4*' sT(1) '+v*(2*' sT(11) '-3*' sT(2) ')+v^2*(2*' sT(12) ' - 2*' ...
    sT(3) ')-' sT(4) '*v^3)*u + v*(' sT(7) '+' sT(8) '*v+' sT(9) '*v^2)+' sT(6)];

end

function strEq = eq2(T)
sT = @(x) ['(' num2str(T(x)) ')'];
strEq = ['(' sT(14) '-' sT(13) '*v)*u^3 + u^2*' '(' sT(11) '+(2*' sT(12) '-2*' sT(10) ...
    ')*v-' sT(11) '*v^2) + u*(' sT(7) '+v*(2*' sT(8) '-3*' sT(6) ')+v^2*(3*' sT(9) ...
    '-2*' sT(7) ') - ' sT(8) '*v^3) + v*(2*' sT(3) '-4*' sT(1) '+v*(3*' sT(4) ...
    '-3*' sT(2) ')+v^2*(4*' sT(5) ' - 2*' sT(3) ')-' sT(4) '*v^3)+' sT(2)];

end

ฉันมีสองคำถาม:

1) เพื่อที่จะใช้ MuPAD ฉันจำเป็นต้องเขียนสมการทั้งสองของฉันใหม่สำหรับระบบสมการเป็นสตริง ดังที่คุณเห็นด้านบน มีวิธีที่ดีกว่าในการทำเช่นนี้ โดยควรไม่มีขั้นตอนสตริงหรือไม่

2) และเกี่ยวกับเอาต์พุตรูปแบบ เมื่อไร

T = [0 0 0 0 0 0 0 0 0 0 1 0 1 0 1];

ผลลัพธ์คือ:

testMupadSolver(T)

ans =

matrix([[u], [v]]) in {matrix([[4.4780323328249527319374854327354], [0.21316518769990291263811232040432]]), matrix([[- 0.31088044854742790561428736573347 - 0.67937835289645431373983117422178*i], [1.1103383836576028262792542770062 + 0.39498445715599777249947213893789*i]]), matrix([[- 0.31088044854742790561428736573347 + 0.67937835289645431373983117422178*i], [1.1103383836576028262792542770062 - 0.39498445715599777249947213893789*i]]), matrix([[0.47897094942962218512261248590261], [-1.26776233072168360314707025141]]), matrix([[-0.83524238515971910583152318717102], [-0.66607962429342496204955062300669]])} union solvelib::VectorImageSet(matrix([[0], [z]]), z, C_)

MuPAD สามารถให้คำตอบเป็นเซตของเวกเตอร์หรือคล้ายกันได้หรือไม่? เพื่อที่จะใช้คำตอบข้างต้น ฉันจำเป็นต้องเรียงลำดับวิธีแก้ปัญหาจากชุดโซลูชันชุดนั้น มีวิธีที่ชาญฉลาดในการทำเช่นนี้หรือไม่? วิธีแก้ปัญหาของฉันจนถึงตอนนี้คือค้นหาสัญญาณที่ฉันรู้ว่าจะมีอยู่ในวิธีแก้ปัญหา เช่น '([[' และเลือกตัวเลขต่อไปนี้ ซึ่งน่าเกลียดมาก และหากวิธีแก้ปัญหาด้วยเหตุผลบางอย่างดูแตกต่างไปจากกรณีที่ฉันทำเล็กน้อย ปกปิดมันใช้งานไม่ได้

แก้ไข

เมื่อฉันใช้โซลูชันที่แนะนำในคำตอบด้านล่างโดย @horchler ฉันจะได้รับโซลูชันเดียวกันกับการใช้งานครั้งก่อน แต่สำหรับบางกรณี (ไม่ใช่ทั้งหมด) อาจใช้เวลานานกว่ามาก เช่น. สำหรับ T ด้านล่าง วิธีแก้ปัญหาที่แนะนำด้านล่างจะใช้เวลามากกว่าหนึ่งนาที ในขณะที่การใช้ evalin (การใช้งานครั้งก่อนของฉัน) จะใช้เวลาหนึ่งวินาที

T = [2.4336 1.4309 0.5471 0.0934 9.5838 -0.1013 -0.2573 2.4830 ...
     36.5464 0.4898 -0.5383 61.5723 1.7637 36.0816 11.8262]

ฟังก์ชั่นใหม่:

function ut1 = testMupadSolver(T)
% # Input T should be a vector of 15 elements
allSolutions = feval(symengine,'numeric::polysysroots', ...
    [eq1(T),eq2(T)],'[u,v]');
end
function eq = eq1(T)
syms u v
eq = -T(13)*u^4 + (4*T(15) - 2*T(10) - T(11)*v)*u^3 + (3*T(13) - 3*T(6) ...
    + v*(3*T(14) -2*T(7)) - T(8)*v^2)*u^2 + (2*T(10) - 4*T(1) + v*(2*T(11) ...
    - 3*T(2)) + v^2*(2*T(12) - 2*T(3)) - T(4)*v^3)*u + v*(T(7) + T(8)*v ...
    + T(9)*v^2) + T(6);
end


function eq = eq2(T)
syms u v
eq = (T(14) - T(13)*v)*u^3 + u^2*(T(11) + (2*T(12) - 2*T(10))*v ...
    - T(11)*v^2) + u*(T(7) + v*(2*T(8) - 3*T(6) ) + v^2*(3*T(9) - 2*T(7)) ...
    - T(8)*v^3) + v*(2*T(3) - 4*T(1) + v*(3*T(4) - 3*T(2)) + v^2*(4*T(5) ...
    - 2*T(3)) - T(4)*v^3) + T(2);
end

มีเหตุผลที่ดีไหมว่าทำไมจึงต้องใช้เวลานานกว่านี้มาก?


person Fija    schedule 30.12.2013    source แหล่งที่มา


คำตอบ (1)


ประการแรก Matlab สื่อสารกับ MuPAD ผ่านคำสั่งสตริง ดังนั้นท้ายที่สุดแล้วจึงไม่มีทางเลี่ยงการใช้สตริงได้ และเนื่องจากเป็นรูปแบบดั้งเดิม หากคุณส่งข้อมูลจำนวนมากไปยัง MuPAD แนวทางที่ดีที่สุดคือการแปลงทุกอย่างเป็นสตริงอย่างรวดเร็วและมีประสิทธิภาพ (sprintf โดยทั่วไปจะดีที่สุด) อย่างไรก็ตาม ในกรณีของคุณ ฉันคิดว่าคุณสามารถใช้ feval แทน evalin ซึ่ง อนุญาตให้คุณส่งผ่านประเภทข้อมูล Matlab ปกติ (ภายใต้ประทุน sym/feval ทำการแปลงสตริงและการโทร evalin) วิธีการนี้จะกล่าวถึงใน บทความ MathWorks สามารถใช้รหัสต่อไปนี้:

T = [0 0 0 0 0 0 0 0 0 0 1 0 1 0 1];
syms u v;
eq1 = -T(13)*u^4 + (4*T(15) - 2*T(10) - T(11)*v)*u^3 + (3*T(13) - 3*T(6) ...
    + v*(3*T(14) -2*T(7)) - T(8)*v^2)*u^2 + (2*T(10) - 4*T(1) + v*(2*T(11) ...
    - 3*T(2)) + v^2*(2*T(12) - 2*T(3)) - T(4)*v^3)*u + v*(T(7) + T(8)*v ...
    + T(9)*v^2) + T(6);
eq2 = (T(14) - T(13)*v)*u^3 + u^2*(T(11) + (2*T(12) - 2*T(10))*v ...
    - T(11)*v^2) + u*(T(7) + v*(2*T(8) - 3*T(6) ) + v^2*(3*T(9) - 2*T(7)) ...
    - T(8)*v^3) + v*(2*T(3) - 4*T(1) + v*(3*T(4) - 3*T(2)) + v^2*(4*T(5) ...
    - 2*T(3)) - T(4)*v^3) + T(2);
allSolutions = feval(symengine, 'numeric::polysysroots',[eq1,eq2],'[u,v]');

อาร์กิวเมนต์สุดท้ายยังต้องเป็นสตริง (หรือละเว้น) และการเพิ่ม ==0 ลงในสมการก็ใช้ไม่ได้ผลเช่นกัน แต่อย่างไรก็ตาม เลขศูนย์ก็ยังมีความหมายโดยนัยอยู่ดี

สำหรับคำถามที่สอง ผลลัพธ์ที่ส่งคืนโดย numeric::polysysroots นั้นไม่สะดวกมาก และไม่ใช่เรื่องง่ายที่จะทำงานด้วย มันคือเซต (DOM_SET) ของเมทริกซ์ ฉันลองใช้ coerce เพื่อแปลงผลลัพธ์เป็นอย่างอื่นโดยไม่เกิดประโยชน์ ฉันคิดว่าคุณดีที่สุดที่จะแปลงเอาต์พุตเป็นสตริง (โดยใช้ char) และแยกวิเคราะห์ผลลัพธ์ ฉันทำสิ่งนี้เพื่อรูปแบบเอาต์พุตที่ง่ายกว่า ฉันไม่แน่ใจว่ามันจะมีประโยชน์หรือไม่ แต่โปรดดูที่ sym2float ซึ่งเพิ่งจัดการเมทริกซ์เชิงสัญลักษณ์ (ส่วน 'matrix([[ ... ]])' ไปที่เอาต์พุตของคุณ) โดยใช้การปรับให้เหมาะสมบางประการ

สิ่งสุดท้าย มีเหตุผลไหมที่ฟังก์ชันตัวช่วยของคุณมีวงเล็บเกินความจำเป็น? ดูเหมือนว่าเพียงพอแล้ว

sT = @(x)num2str(T(x),17);

or

sT = @(x)sprintf('%.17g',T(x));

โปรดทราบว่า num2str จะแปลงเป็นทศนิยมสี่ตำแหน่งตามค่าเริ่มต้นเท่านั้น int2str (หรือ %d ควรใช้หาก T(x) เป็นจำนวนเต็มเสมอ)

person horchler    schedule 31.12.2013
comment
สวัสดีและขอบคุณ! ฉันจะเริ่มต้นด้วยการตอบความคิดเห็นล่าสุดของคุณ วงเล็บที่ไม่จำเป็นนั้นมีความจำเป็นเนื่องจากองค์ประกอบใน T ไม่จำเป็นต้องเป็นค่าบวก ดังนั้นในเวอร์ชันปัจจุบันของฉันจึงจำเป็น แต่บางทีอาจไม่ตรงกับคำแนะนำของคุณ สำหรับข้อเสนอแนะอื่นๆ ของคุณ พวกเขาดูดีจริงๆ แต่ฉันไม่สามารถทดสอบได้อย่างถูกต้องจนกว่าฉันจะกลับมาทำงานภายในสองสามวัน - person Fija; 31.12.2013
comment
ปล. ฉันมีชื่อเสียงน้อยเกินไปที่จะโหวตคำตอบของคุณ แต่จะทำเช่นนั้นโดยเร็วที่สุด ดีเอส - person Fija; 31.12.2013
comment
โปรดดูการแก้ไขของฉัน @horchler ฉันได้ปฏิบัติตามคำแนะนำของคุณแล้วและใช้งานได้ดี ยกเว้นบางครั้งอาจใช้เวลานานกว่านั้นมาก - person Fija; 03.01.2014
comment
@Fija: คุณแน่ใจหรือว่า evalin ใช้เวลาไม่นานเท่ากัน? กล่องเครื่องมือสัญลักษณ์จะแคชผลลัพธ์ ดังนั้นหากคุณรัน feval ในสมการใหม่ก่อน คุณอาจไม่ได้รับการเปรียบเทียบที่ยุติธรรม คุณสามารถลองโทร clear all และ/หรือ clear classes อาจมีเรื่องอื่นเกิดขึ้น ไม่ว่าในกรณีใด คุณสามารถคาดหวังได้ว่า feval จะช้าลงเล็กน้อยด้วยเหตุผลที่ระบุไว้ในคำตอบของฉัน - person horchler; 03.01.2014
comment
ไม่ ฉันแน่ใจว่าทำทั้ง clear all และ clear classes ระหว่างการวิ่งและ feval ใช้เวลาประมาณ 70 วินาที และ evalin ประมาณ 2 วินาที แต่ฉันคิดว่าคำตอบของคุณดีและนี่เป็นปัญหาอื่น ดังนั้นฉันจะยอมรับคำตอบของคุณ =) - person Fija; 04.01.2014
comment
@Fija: สิ่งที่อาจเกิดขึ้นคือในบางกรณีโซลูชันการวิเคราะห์ที่ซับซ้อนอาจถูกส่งคืนในขั้นตอนกลางซึ่งมีราคาแพงในการประเมินเชิงตัวเลข เป็นไปได้ว่าวิธีที่คุณระบุรหัส evalin จะทำให้โซลูชันตัวเลข (คลาสสัญลักษณ์ แต่เป็นตัวเลขทั้งหมด) ที่จะได้รับการประเมินโดยตรง อย่างหลังอาจเร็วกว่ามากหากสูตรที่ส่งคืนโดย feval มีราคาแพง คุณอาจดูว่าคุณสามารถใช้ vpa กับ T หรือสมการของคุณได้หรือไม่ - person horchler; 04.01.2014
comment
ฉันคิดว่ามันอาจจะตรงกันข้ามเช่นกัน และในบางกรณี วิธีแก้ปัญหาเชิงตัวเลขได้รับการประเมินซึ่งสุดท้ายก็ช้า ประเด็นก็คือสตริงที่ feval สร้างมีแนวโน้มว่าจะแตกต่างในแง่ของค่าคงที่/พารามิเตอร์ตัวเลข การรวม/แยกจุดทศนิยม (.0) สามารถเปลี่ยนวิธีประเมินสิ่งต่างๆ ภายในได้ - person horchler; 04.01.2014