1樓:化學工程
clc;clear
t=1790:10:2000;
x=[3.9 5.3 7.
2 9.6 12.9 17.
1 23.2 31.4 38.
6 50.2 62.9 76.
0 92.0 106.5 123.
2 131.7 150.7 179.
3 204.0 226.5 251.
4 281.4];
y=log(x);
p=polyfit(t,y,1)
r=p(1),x0=exp(p(2))
y=polyval(p,t);
x=exp(y);
%方法二,非
線性迴歸
fun=inline('a(1)*exp(a(2)*t)','a','t')
format short g
a=nlinfit(t,x,fun,[1.7139e-010 0.02]);
x0=a(1),rr=a(2)
x1=a(1)*exp(a(2)*t);
plot(t,x,'o',t,x,t,x1)結果:p =
0.020219 -34.393r =0.020219
x0 =
1.1565e-015
x0 =
1.3155e-010
rr =
0.014223
從圖形中可以看版
出,非線性方法效權果更好。
matlab求解微分方程並畫圖 20
2樓:李修靈
看標題以為你要求微分方程吶,結果是畫dr/dr vs. r
% 畫出圖中的公式
% 定義微分方程函式
drdr = @(r) 0.89 ./ r .* exp(-(log(r) + 0.84).^2 / 0.086);
% 在(0, 10]上畫圖
r = linspace(0.01, 3, 500);
dr = drdr(r);
figure(10);
plot(r, dr, '^b', 'marke***cecolor', 'blue', 'displayname', 'dr/dr vs r');
xlabel('r', 'fontsize', 16);
ylabel('dr/dr', 'fontsize', 16);
legend('show');
作圖的結果是
這是一個一階線性微分方程, 可以用龍格庫塔求解. 首先用一個.m檔案
定義圖中的微分方程. 然後再用matlab的ode45函式求解這個方程.
舉個例子, 建立一個solvef**.m的檔案如下
% 呼叫matlab的`ode45`函式實現求解. 具體說明在
% matlab命令列中輸入`doc ode45`查詢.
% 簡單來說, 需要兩步. 1) 微分方程定義; 2) ode45引數設定和呼叫.
% 具體如下.
% 設定ode45引數, 並呼叫求解微分方程`dr = odefunc(r, r)`. 該方程定義
% 見後一個函式.
function solution = solvef**(initial_value)
% initial_value = [r0, r0];
r0 = initial_value(1);
r0 = initial_value(2);
% 由於提供的初值[r0 r0] = [0.43 0.5]; 或 [0.55 0.9]. 如果求解區間從(0, 10]
% 那麼需要分成(0 r], [r, 10]兩部分求解.
% [r0, 10]上求解
sol1 = ode45(@odefunc, [r0 10], r0);
% (0, r0]上求解
sol2 = ode45(@odefunc, [r0 0.01], r0); % log(r)在0處inf, 所以用0.01代替.
% sol結構體是微分方程的解, 包含`r` - 自變數, `r` - 要求解的函式r(r)
% 的數值. 其它r對應的r值, 可以用 `r = deval(sol, r);`計算. 其結果與
% sol中儲存的r, r(r)相對精度一致.
solution.sol1 = sol1;
solution.sol2 = sol2;
end% 定義微分方程. 用dr表示dr/dr. 一階線性微分方程的一般形式為dr = f(r, r)
% 儘管題目微分方程右邊未出現r, 此變數也不可省略.
function dr = odefunc(r, r)
dr = 0.89 ./ r .* exp(-(log(r) + 0.84).^2 / 0.086);
end這個函式的作用就是求解最後幾行定義的 dr/dr = f(r, r)這個微分方程. 解存在sol結構體中.(這個例子是在sol1和sol2裡面)
怎麼用這個sol解具體算每個r處的r呢, 需要用到deval. 可以再建立一個main_0.m檔案, 內容為
% 呼叫函式sovlef**. 求解微分方程
% 初值 r= 0.43, r = 0.5
s1 = solvef**( [0.43 0.5] );
% 繪圖
r1 = linspace(0.01, 0.43, 50);
r1 = deval(s1.sol2, r1);
r2 = linspace(0.43, 2, 50);
r2 = deval(s1.sol1, r2);
r = [r1, r2]; % 在 0.01 - 2上的解
r = [r1, r2]; % 在 0.01 - 2上的解
figure(1); hold on;
plot(r, r, 'b-', 'linewidth',2, 'displayname', 'r=0.43, r=0.5');
% 初值 r=0.55, r = 0.9
s2 = solvef**( [0.55 0.9] );
% 繪圖
r1 = linspace(0.01, 0.55, 50);
r1 = deval(s2.sol2, r1);
r2 = linspace(0.55, 2, 50);
r2 = deval(s2.sol1, r2);
r = [r1, r2]; % 在 0.01 - 2上的解
r = [r1, r2]; % 在 0.01 - 2上的解
figure(1); hold on;
plot(r, r, 'r-', 'linewidth',2, 'displayname', 'r=0.55, r=0.9');
legend('show')
xlabel('r', 'fontsize', 16);
ylabel('\int_0^r dr', 'fontsize', 16)
那麼這個微分方程在不同的初值條件下的解, 就可以畫成下圖
3樓:華凡劍春竹
如何用matlab求解微分方程並畫圖,可以先用dsolve()或ode()求出其微分方程(組)的解析解或數值解,然後用plot()繪製其圖形。
例如:解微分方程 y'=y-2t/y,y(0)=1,0 3、當然嘍,使用dsolve()或ode()求解要根據題意去分析,來決定用那個函式。一般來說,用ode45求解微分方程(組)的數值解用點比較多。 4樓:莫悟軒轅良俊 syms tv=dsolve('dv=(190.708-90.64*v^2)/47.27','v(0)=0','t'); t=0:0.00001:0.002; v=eval(v); plot(t,v) 使用這樣的方法求解,但從結果看好像你的方程有問題! 如果說是範圍的話應該是滿足與三角形三條邊直線方程有關的不等式組但是如果你有三版點a x1,y1 b x2,y2 c x3,y3 組成三角形權想知道某點 x,y 是否在三角型裡面可以用matlab函式 in on inpolygon x,y,x1 x2 x3 y1 y2 y3 返回的in和on都是邏輯... n0 1 n1 200 取最左邊的一個週期的邊界t0 200 週期for i 0 4 5個週期t n0 0.1 n1 y 0.002 1.0191 t t0 i 表示式 plot t,y hold on n0 n0 t0 左邊界右移一個週期n1 n1 t0 右邊界右移一個週期end 第6個週期 t ... subplot就是將figure中的影象劃分為幾塊，每塊當中顯示各自的影象，有利於進行比較。比如example裡面有這樣的例子 in e 3.2 4.1 5.0 5.6 outgo 2.5 4.0 3.35 4.9 subplot 2,1,1 plot in e subplot 2,1,2 plot...matlab如何用inpolygon函式判斷點是否在園內
如何用MATLAB畫周期函式，如何用MATLAB畫周期函式？
如何用matlab中subplot的使用