Menemukan semua solusi sistem persamaan non-linier dengan MuPAD

Pertanyaan saya adalah apakah ada cara yang baik untuk menggunakan fungsi MuPAD dalam skrip Matlab. Latar belakangnya adalah saya mempunyai masalah dimana saya harus mencari semua solusi dari himpunan persamaan non-linier. Solusi sebelumnya adalah dengan menggunakan solve di Matlab, yang berfungsi untuk beberapa simulasi saya (yaitu, beberapa set input T) tetapi tidak selalu. Jadi saya menggunakan MuPAD dengan cara berikut:

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

Saya punya dua pertanyaan:

1) Untuk menggunakan MuPAD saya perlu menulis ulang dua persamaan saya untuk sistem persamaan sebagai string, seperti yang Anda lihat di atas. Apakah ada cara yang lebih baik untuk melakukan ini, sebaiknya tanpa langkah string?

2) Dan mengenai format keluaran; Kapan

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

hasilnya adalah:

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_)

Bisakah MuPAD memberikan solusi sebagai himpunan vektor atau sejenisnya? Untuk menggunakan jawaban di atas, saya perlu memilah solusi dari rangkaian solusi tersebut. Apakah ada cara cerdas untuk melakukan ini? Solusi saya sejauh ini adalah menemukan tanda-tanda yang saya tahu akan ada dalam solusi, seperti '([[' dan memilih angka-angka berikut, yang benar-benar jelek, dan jika solusi karena alasan tertentu terlihat sedikit berbeda dari kasus yang saya miliki menutupinya tidak berhasil.

Sunting

Ketika saya menggunakan solusi yang disarankan dalam jawaban di bawah ini oleh @horchler, saya mendapatkan solusi yang sama seperti implementasi saya sebelumnya. Namun untuk beberapa kasus (tidak semua) memerlukan waktu lebih lama. Misalnya. untuk T di bawah solusi yang disarankan di bawah ini membutuhkan waktu lebih dari satu menit sementara menggunakan evalin (implementasi saya sebelumnya) membutuhkan waktu satu detik.

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]

Fungsi baru:

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

Apakah ada alasan bagus mengapa dibutuhkan waktu lebih lama?


matlab symbolic-math mupad
person Fija    schedule 30.12.2013    source sumber

Jawaban (1)


arrow_upward
2
arrow_downward

Pertama, Matlab berkomunikasi dengan MuPAD melalui perintah string sehingga pada akhirnya tidak ada cara untuk menyiasati penggunaan string. Dan karena ini adalah format asli, jika Anda meneruskan data dalam jumlah besar ke MuPAD, pendekatan terbaik adalah mengonversi semuanya menjadi string dengan cepat dan efisien (sprintf biasanya yang terbaik). Namun, dalam kasus Anda, menurut saya Anda dapat menggunakan feval alih-alih evalin yang memungkinkan Anda meneruskan tipe data Matlab biasa (di balik terpal sym/feval melakukan konversi string dan memanggil evalin). Metode ini dibahas dalam artikel MathWorks. Kode berikut dapat digunakan:

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]');

Argumen terakhir masih harus berupa string (atau dihilangkan) dan menambahkan ==0 ke persamaan juga tidak berfungsi, tetapi angka nol tetap tersirat.

Untuk pertanyaan kedua, hasil yang dikembalikan oleh numeric::polysysroots sangat merepotkan dan tidak mudah untuk dikerjakan. Ini adalah himpunan (DOM_SET) matriks. Saya mencoba menggunakan coerce untuk mengonversi hasilnya ke hasil lain, tetapi tidak berhasil. Saya pikir Anda sebaiknya mengonversi output menjadi string (menggunakan char) dan mengurai hasilnya. Saya melakukan ini untuk format keluaran yang lebih sederhana. Saya tidak yakin apakah ini akan membantu, tapi jangan ragu untuk melihat sym2float yang hanya menangani matriks simbolis (bagian 'matrix([[ ... ]])' adalah keluaran Anda) menggunakan beberapa pengoptimalan.

Hal terakhir. Apakah ada alasan mengapa fungsi pembantu Anda menyertakan tanda kurung yang berlebihan? Ini tampaknya cukup

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

or

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

Perhatikan bahwa num2str hanya mengonversi ke empat tempat desimal secara default. int2str (atau %d harus digunakan jika T(x) selalu berupa bilangan bulat).

person horchler    schedule 31.12.2013
comment
Hai, dan terima kasih! Saya akan mulai dengan menjawab komentar terakhir Anda, tanda kurung yang berlebihan diperlukan karena elemen di T belum tentu positif sehingga dalam versi saya saat ini elemen tersebut diperlukan, tetapi mungkin tidak sesuai dengan saran Anda. Mengenai saran Anda yang lain, saran tersebut terlihat sangat bagus, tetapi saya tidak dapat mengujinya dengan benar sampai saya kembali bekerja dalam beberapa hari. - person Fija; 31.12.2013
comment
hal. Saya memiliki reputasi yang terlalu kecil untuk memberi suara positif pada jawaban Anda, tetapi akan melakukannya sesegera mungkin. ds. - person Fija; 31.12.2013
comment
Silakan lihat hasil edit saya @horchler. Saya telah menerapkan saran Anda dan berfungsi dengan baik, hanya saja terkadang memerlukan waktu lebih lama. - person Fija; 03.01.2014
comment
@Fija: Apakah Anda yakin evalin tidak memakan waktu lama? Kotak alat simbolis menyimpan hasil dalam cache, jadi jika Anda menjalankan feval pada persamaan baru terlebih dahulu, Anda mungkin tidak mendapatkan perbandingan yang adil. Anda dapat mencoba menelepon clear all dan/atau clear classes. Mungkin ada hal lain yang terjadi. Bagaimanapun, Anda dapat mengharapkan feval menjadi sedikit lebih lambat karena alasan yang diuraikan dalam jawaban saya. - person horchler; 03.01.2014
comment
Tidak, saya yakin, melakukan clear all dan clear classes antara proses dan feval membutuhkan waktu sekitar 70 detik dan evalin sekitar 2 detik. Tapi menurut saya jawaban Anda bagus dan ini masalah yang berbeda jadi saya akan menerima jawaban Anda. =) - person Fija; 04.01.2014
comment
@Fija: Apa yang mungkin terjadi adalah dalam beberapa kasus solusi analitik yang rumit dapat dikembalikan pada langkah perantara yang mahal untuk dievaluasi secara numerik. Mungkin saja cara Anda menentukan kasus kode evalin solusi numerik (kelas simbolik tetapi seluruhnya numerik) untuk dievaluasi secara langsung. Yang terakhir bisa lebih cepat jika rumus yang dikembalikan oleh feval mahal. Anda mungkin melihat apakah Anda dapat menggunakan vpa pada T atau persamaan Anda. - person horchler; 04.01.2014
comment
Saya kira bisa juga sebaliknya dan dalam beberapa kasus solusi numerik dievaluasi namun akhirnya menjadi lambat. Intinya adalah string yang feval buat kemungkinan besar berbeda dalam hal konstanta/parameter numerik. Penyertaan/pengecualian koma desimal (0,0) dapat mengubah cara penilaian internal. - person horchler; 04.01.2014