1. 程式人生 > 其它 >《白色電路2》——RLC一階與二階電路的冬季戀歌

《白色電路2》——RLC一階與二階電路的冬季戀歌

技術標籤:matlab電磁學

寒冷的風中傳來MATLAB報錯的聲音,它讓我們原本凌亂的三種旋律緊緊聯絡到了一起。公式法,deslove,ode,三種可以實現的方法,RC,RL,RLC多種不同的電路,雪菜,冬馬,小春,無論是誰都付出了全心全意。無論是誰都在以堅強意志向前。
實現思路:
對於RC,RL,RLC電路而言,自動實現的程式關鍵在於列出並解決微分方程,我們既可以直接選用課本上的公式代入計算,也可以用deslove函式求出解析解,用ode方法也可以求出數值解,北原春希又會選擇誰呢?
(1)、使用解析解desolve方法
選擇了desolve方法就如同選擇了最為溫柔可愛的小木曾雪菜,她或許不是那麼的全能,但是卻可以讓你很簡單地解決很多問題,在解決此次問題時,我首先選用的即使這種方法,在寫入app designer前,我首先在MATLAB裡寫出基本程式碼,然後再轉入app designer中。

Us=input('please input the Us:');
        R=input('please input the R:');
        C=input('please input the C:');
        L=input('please input the L:');
        Uc0=input('please input the Uc:');
        Il0=input('please input the Il:');
        syms uc(t)
        eqn = L*C*diff(uc,t,2)+  R*C*diff(
uc,t) + uc == Us; Duc=diff(uc,t); cond=[uc(0)==Uc0,Duc(0)==-Il0/C]; Uc(t) = dsolve(eqn,cond) t=0:0.0001:0.1; Uc1=Uc(t); plot(t,Uc1);

但是雪菜的溫柔可愛不代表她沒有心機,不懂事。在我將此程式碼寫至app designer時,我發現在app designer中使用syms有時似乎會遇見問題,因此採用了更為傳統的desolve的用法類似下例:

y=desolve('Dy=1')

對於字母使用deval即可。

但是溫柔可愛的雪菜似乎從來不會進入春希同學的心中,即使使用了更為傳統的desolve用法,我發現在一些情況下會出現結果為NaN的無解情況。是的,雪菜或許很溫柔,很可愛,但是,有些事情,她真的無能為力,她可以在春希同學因為苦苦等待和紗同學而感冒生病時照顧他,可以一直原諒他,但是這場愛戀如今依然寸步難行。此時,我選擇了另一種方法。
(2)、使用ode求數值解。
春希同學或許永遠不會知道,和紗同學看似冷峻的外表下的動搖與希冀,就像他一直以來不知道那個在隔壁用鋼琴聲鼓勵他的人就是和紗同學一樣。在用desolve面對一些情況出現NaN的無解之後,我選擇了ode方法,這種方法雖然可能在理解上較為困難,但是卻保證有解。就像和紗同學,雖然可能外表冷峻,但是她卻最適合春希。
在使用ode時程式碼如下:

case 'RLC'
        R1=app.R.Value;
        C1=app.C.Value;
        L1=app.L.Value;
        Us=app.S.Value;
        Uc0=app.Z1.Value;
        Il0=app.Z2.Value;
        t1=20*L1/R1;
        app.T21=t1;
        tspan=[0 t1];
        app.sol = ode45(@(t,y)[y(2); (Us-y(1)-R1*C1*y(2))/(L1*C1)] ,tspan,[Uc0; Il0]);
        x = linspace(0,t1,250);
        y = deval(app.sol,x);
        app.y1 = -L1*C1*diff(y(2,:))/(t1/250);
        app.y1=[app.y1 app.y1(249)];
        plot(app.UIAxes,x,y(1,:),'b',x,-y(2,:)*R1*C1,'r',x,app.y1,'g');
        title(app.UIAxes,'RLC串聯電路');
        ylabel(app.UIAxes,'U/V');
        legend(app.UIAxes,'uc','ur','ul');
        name='tm';
        namet=t1;
        k={ '\' name namet};
        app.UITable.Data=k;

(3)、公式法
步入大學後,和紗前往歐洲,春希與雪菜貌合神離,如果面對剪不斷理還亂的困境,不如抽身而出。偶然認識的小春,也很可愛不是嗎?在用數值解和解析解做出答案後,我想嘗試跳出微分方程求解的圈子,直接代入公式求解,就像春希未曾與小春度過高中三年一樣,這一過程必定充滿困難與摩擦。
首先第一點在於能夠完整求出各種情況下的公式,書上所給公式中均默認了一階倒對應的初值為0,因此我們需要從頭推出一遍情況下uc,il均不為零情況下的公式。
該處求微分方程過程過於複雜略去。
結果為推理和驗證了書上該章所有公式,並且對書上未給出的公式進行了討論,下面代入計算即可。

 case 'RLC'
        R1=app.R.Value;
        C1=app.C.Value;
        L1=app.L.Value;
        Us=app.S.Value;
        Uc0=app.Z1.Value;
        Il0=app.Z2.Value;
        if R1>2*(L1/C1)^(1/2)
            p1=-R1/(2*L1)+((R1/(2*L1))^2-1/(L1*C1))^(1/2);
            p2=-R1/(2*L1)-((R1/(2*L1))^2-1/(L1*C1))^(1/2);
            A1=(p2*(Uc0-Us)*C1+Il0)/((p2-p1)*C1);
            A2=(p1*(Uc0-Us)*C1+Il0)/((p1-p2)*C1);
            t1=log(p2/p1)/(p1-p2);
            app.T21=t1;
            t=0:t1/1000:10*t1;
            app.uc=A1*exp(p1*t)+A2*exp(p2*t)+Us;
            app.il=-C1*diff(app.uc)*(1000/t1);app.il=[app.il(1) app.il];
            app.ur=app.il*R1;
            app.ul=L1*diff(app.il)*(1000/t1);app.ul=[app.ul(1) app.ul];           
        end
        if R1==2*(L1/C1)^(1/2)
            p=-R1/(2*L1);
            t1=2*L1/R1;
            app.T21=t1;
            t=0:t1/1000:10*t1;
            app.uc=(Uc0-Us+(-Il0/C1-p*(Uc0-Us))*t)*exp(p*t)+Us;
            app.il=-C1*diff(app.uc)*(1000/t1);app.il=[app.il(1) app.il];
            app.ur=app.il*R1;
            app.ul=L1*diff(app.il)*(1000/t1);app.ul=[app.ul(1) app.ul];
        end
        if R1<2*(L1/C1)^(1/2)
            w=(1/(L1*C1)-(R1/(2*L1))^2)^(1/2);
            q=R1/(2*L1);
            b=atan(w/q);
            p1=-R1/(2*L1)+((R1/(2*L1))^2-1/(L1*C1))^(1/2);
            p2=-R1/(2*L1)-((R1/(2*L1))^2-1/(L1*C1))^(1/2);
            t1=abs(log(p2/p1)/(p1-p2));
            app.T21=t1;
            t=0:t1/1000:10*t1;
            app.uc=(Uc0-Us).*(q^2+w^2)^(1/2).*exp(-q.*t).*sin(w.*t+b)/w-Il0.*exp(-q.*t).*sin(w.*t)./(w.*C1)+Us;
            app.il=-C1*diff(app.uc)*(1000/t1);app.il=[app.il(1) app.il];
            app.ur=app.il*R1;
            app.ul=L1*diff(app.il)*(1000/t1);app.ul=[app.ul(1) app.ul];
        end
        plot(app.UIAxes,t,app.uc,'b',t,app.ur,'r',t,app.ul,'g');
        title(app.UIAxes,'RLC串聯電路');
        ylabel(app.UIAxes,'U/V');
        legend(app.UIAxes,'uc','ur','ul');
        name='tm';
        namet=t1;
        k={ '\' name namet};
        app.UITable.Data=k;

當然,可愛的小春也會有很多困擾人的地方,例如在∆<0時對複數的處理等,在此處我直接進行取模處理。
(4)、結果
《白色相簿2》的目的不在於讓你抉擇,無論你是選擇了雪菜,冬馬還是小春你都會幸福,就像在我使用完這三種方法後,我對這三種方法求出的值進行了比較,雖然誤差存在,但是相對較小,幸福有很多種,但是卻一定會幸福。
圖中為用ode和公式法做出的結果對比。圖中為用ode和公式法做出的結果對比。
(5)、各種方法的優劣對比
在這裡插入圖片描述(6)、關於圖中總時間的選擇
RC,RL電路中:選用4倍時間常數,此時可以較好地反應變化情況;
RLC,GLC電路中:選用震盪達到最大值時間的10倍,此時可較好的反應震盪和變化情況。
(7)、關於數值傳遞方法的選擇
在本例中,傳遞數值的方法有兩種,一種為傳遞函式關係,一種傳遞值對應矩陣,傳遞函式關係需要重新進行計算,但是對空間佔用小,傳遞值矩陣查詢效果更快,但是對空間佔用大,綜合比較兩方法後,我們發現該兩種方法相對誤差控制在 1%以下。
在求解過程中遇到的問題:
(1)、傳達不了的戀情,已經不需要了。因為已經不再有人,值得去愛了。
在剛開始時一直困擾我的一點在於在不同回撥中傳遞變數,但是後來發現可以直接在開頭定義全域性變數即可。。。
(2)、這兩份快樂交織在一起,我應該更快樂才對,但是。。。
Desolve現在似乎鼓勵先用syms然後再求數值解,但是這種被鼓勵的方法在明明應該很高效的app designer中可能會失效。
(3)、對不起,和紗…雖然我借的時間有些長,春希君,還給你吧。
還有另一種方法在不同回撥間傳遞變數,可以將變數暫時放在UItable中,既能夠顯示,還能夠暫時儲存已讓其他回撥使用。
(4)、為什麼,同意我呆在你身邊;為什麼,不能接受我啊
UItable中的格式為cell直接用()讀取的是cell,只有用{}才能讀取到UItable中的值。
(5)、為什麼你會這麼熟練啊!你到底要把我甩開多遠你才甘心啊!?
App designer中不會把變數值放在工作區,對bug的查詢產生困難。