(هادی)
کاربر فعال تالار ریاضی ,
احتمالا خطای عددی باشه. با عوض کردن دستور حل معادله دیفرانسیل، یا با کوچک کردن گام زمانی حل، باید درست بشهمشاهده پیوست 258620
مثل این...در واقع باید آخرش steady بشه ولی نوسانی در میاد
احتمالا خطای عددی باشه. با عوض کردن دستور حل معادله دیفرانسیل، یا با کوچک کردن گام زمانی حل، باید درست بشهمشاهده پیوست 258620
مثل این...در واقع باید آخرش steady بشه ولی نوسانی در میاد
گاهی اوقات روشهای ضمنی چنین جوابی میدن. اگه روشهای صریح به جواب همگرا بشن، جوابشون احتمالاً چنین مشکلی نخواهد داشت. واسه همین، اگه ode15s استفاده کردید، به جاش، از ode45 و ode23 رو تست کنید ببینید جواب همگرا میشه یا نه.
از روش ضمنی استفاده کردم و آخرt,c) plot) را نوشتم.از ode استفاده نکردم
احتمالا خطای عددی باشه. با عوض کردن دستور حل معادله دیفرانسیل، یا با کوچک کردن گام زمانی حل، باید درست بشه
توی اون پست قبلی هم گفتم از روش صریح استفاده کنید.
اگه می خواید حتماً از روش ضمنی استفاده کنید، چند تا کار که در CFD مرسومه می تونید انجام بدید:
1. با استفاده از تحلیل پایداری (مثلاً روش نیومن) بازه یا مقدار مناسب گامها رو پیدا کنید. پیشنهاد من اینه که به جای این کار برنامه رو با جایگشتهای مختلف گامهای dr، dz و dt اجرا کنید و هر جوابی که بهتر بود رو نگه دارید. توجه کنید که در اینگونه مسائل کوچک کردن یک گام لزوماً منجر به جواب بهتری نمی شه؛ ممکنه بزرگتر کردن اون گام منجر به جواب بهتری بشه.
2. اگه همه این کارها رو انجام دادید و نتیجه مطلوب رو به دست نیاوردید، می تونید جواب رو از یک فیلتر low pass عبور بدید؛ البته اگه این کار در حین حل معادلات انجام بشه بهتره تا بعد از حل معادلات.
با سلام .
نماز و روزه و طاعات و عبادات شما عزیزان قبول درگاه حق انشاالله.
میخواستم این کد رو برام تصحیح بکنید که بتونم اعداد تو جدول رو بدست بیارم،
یه دستگاه دو معادله ای هستش که باید با روش رانگ کوتا مرتبه 4 حل بشه که دو مدل داره یکی اینکه رانگ کوتا رو به صورت یه تابع تعریف بکنیم و نو برنامه ای که مینویسیم فقط معادلات رو جایگذرای میکنیم و یکی دیگه اینکه تو برنامه اصلی معادلات رو مینویسیم و تو برنامه ای که داریم تمام k ها رو حساب میکنیم.
من الان این مدل دومی رو نیاز دارم.
نمیخوام پیشرفته باشه، مثل همین که خودم نوشتم فقط تصحیح بشه .
با تشکر فراوان.
http://s7.picofile.com/file/8257314892/kode_2_moadeleie.rar.html
سلام
نماز و روزه شما هم قبول باشه
به یک همچین فرمی باید بنویسید:
clc
h = 0.5;
x0 = 0;
y = [4; 6];
for i=1:4;
x=x0+i*h;
k1=RK(t,y);
% k2=...;
% k3=...;
% k4=...;
y=y+(h/6)*(k1+2*k2+2*k3+k4)
end
و ورودی تابع دینامیک هم باید درست بشه:
function eq=RK(x, y);
eq = [-0.5*y(1);4-0.3*y(2)-0.1*y(1)];
موفق باشید
پس دو تا تابع f1 و f2 هم نیاز داری که در واقع از هم جدا بشن
مهندس یعنی باید 2 تا function درست بکنم و داخل هر کدوم جداگونه تابع ها رو تعریف بکنم؟
باز برمیگردم به خونه اول فکر کنم!
من همون روز که بهمون گفت چی میخواد از ما 2 تا function نوشتم و یه فایل اسکریپت و وقتی میخواستم 8 تا k ها رو حساب بکنم برای هر 4 تاش به یکیشون لینک میدادم ولی استاد گفت نه این خیلی مبتدیانه هستش و همه اش باید تو یک function باشه!
................
یا اینکه منظورتون اینه که این 2 تا f1 و f2 رو تو خود فایل اسکریپتی که دارم مینویسم تعریف بکنم و نه داخل function ؟
والله به خدا به عقل من هم جور در نمیاد !بالاخره یا باید متغیرها رو از هم جدا کنی و هر کدوم رو جداگانه پهاربار صدا بزنی، یا همه رو یک کاسه کنی توی یک تابع
اگه همه رو یکی کردی (که منطقی تر هم همینه)، همون تابع رو چهار بار صدا میزنی ...
به عقلم جور در نمیاد که همه رو بریزی توی یک تابع، ولی به جای چهار بار، هشت بار صدا بزنی اون تابع رو
dt .dz.dr را تغییر دادم.هرچی کمتر باشن جوابا بهترن...ببخشید low pass چی هست؟چجوری بکار میره؟
>> help sgolay
خب شما اگه f1 در حکم ydot هست و f2 در حکم tdot، ...
باید k1 و k2 و k3 و k4 رو فقط با y جمع کنید و اون چهارتای دیگه رو فقط با t
منظورم اینه کهt+0.5*h*k1 مشکل داره
با سلام .
نماز و روزه و طاعات و عبادات شما عزیزان قبول درگاه حق انشاالله.
میخواستم این کد رو برام تصحیح بکنید که بتونم اعداد تو جدول رو بدست بیارم،
یه دستگاه دو معادله ای هستش که باید با روش رانگ کوتا مرتبه 4 حل بشه که دو مدل داره یکی اینکه رانگ کوتا رو به صورت یه تابع تعریف بکنیم و نو برنامه ای که مینویسیم فقط معادلات رو جایگذرای میکنیم و یکی دیگه اینکه تو برنامه اصلی معادلات رو مینویسیم و تو برنامه ای که داریم تمام k ها رو حساب میکنیم.
من الان این مدل دومی رو نیاز دارم.
نمیخوام پیشرفته باشه، مثل همین که خودم نوشتم فقط تصحیح بشه .
با تشکر فراوان.
http://s7.picofile.com/file/8257314892/kode_2_moadeleie.rar.html
function f = RK(x,y)
f = [-0.5*y(1); 4-0.3*y(2)-0.1*y(1)];
% Runge8
clear('all'), close('all'), clc
h=0.5;
x0=0;
xf = 2;
x = x0 : h : xf; x = x(:);
y0 = [4, 6];
y = zeros(length(x), length(y0));
y(1,:) = y0(:).';
for i=1:length(x)-1
k1=RK(x(i),y(i,:));
ytmp = y(i,:) + k1.'*h/2;
k2=RK(x(i)+h/2,ytmp);
ytmp = y(i,:) + k2.'*h/2;
k3=RK(x(i)+h/2,ytmp);
ytmp = y(i,:) + k3.'*h;
k4=RK(x(i)+h,ytmp);
y(i+1,:)=y(i,:)+(h/6)*(k1+2*k2+2*k3+k4).';
end
xy = [x, y]
plot(x,y,'o-')
در مورد روزه؛ از 12 ماه سال این تنها ماهی هست که میشه درش روزه خواری کرد، من حیفم میاد چنین فرصتی رو از دست بدم.
اما جواب سؤال شما:
کد:function f = RK(x,y) f = [-0.5*y(1); 4-0.3*y(2)-0.1*y(1)];
کد:% Runge8 clear('all'), close('all'), clc h=0.5; x0=0; xf = 2; x = x0 : h : xf; x = x(:); y0 = [4, 6]; y = zeros(length(x), length(y0)); y(1,:) = y0(:).'; for i=1:length(x)-1 k1=RK(x(i),y(i,:)); ytmp = y(i,:) + k1.'*h/2; k2=RK(x(i)+h/2,ytmp); ytmp = y(i,:) + k2.'*h/2; k3=RK(x(i)+h/2,ytmp); ytmp = y(i,:) + k3.'*h; k4=RK(x(i)+h,ytmp); y(i+1,:)=y(i,:)+(h/6)*(k1+2*k2+2*k3+k4).'; end xy = [x, y] plot(x,y,'o-')
مهندس جان یه دنیا ممنون،
بسیار تر و تمیز و پیشرفته!
البته اگر کمی مبتدیانه تر هم میشد حل بشه عالی بود.
فقط یه مورد بنده دوست دارم مدلی که چندین k رو استفاده میکنیم رو یادبگیرم!چونکه خیلی ضروریه.
اخه قراره تو امتحان ده معادله ای بهمون بده که میگه 40 تا k باید براش تعریف بکنیم!
لطفا اگر زحمتی نیست میشه این رو به صورت 8 تا k برای y1 ,y2 کد نویسی بکنید ؟
سپاسگذار شما میشیم.
k1=RK(x(i),y(i,:));
k11 = k1(1);
k12 = k1(2);
ytmp1 = y(i,1) + k11*h/2;
ytmp2 = y(i,2) + k12*h/2;
اون کار دیگه با توجه به تعریف جمع ماتریسها و ضرب اسکالر در یک ماتریس کار ساده ای هستش. مثلاً برای k1 و ytmp می تونید بنویسید:
بقیه هم مشابه به همینه.
کد:k1=RK(x(i),y(i,:)); k11 = k1(1); k12 = k1(2); ytmp1 = y(i,1) + k11*h/2; ytmp2 = y(i,2) + k12*h/2;
می تونید تابع RK رو هم تغییر بدید تا به جای اینکه در خروجی یه بردار دو درایه ای بده، دو تا اسکالر بده؛ این طوری شباهتش به برنامه های متلب کمتر میشه.
بعدش هم برای امتحان این روش که معقول نیست؛ سر امتحان باید یه ماشین حساب مهندسی مثلاً CASIO fx-4500P همراهتون داشته باشید و برای این جور محاسبات برنامه بنویسید.
ممنون مهندس جان ، یه خورده سخته ، دارم روش کار میکنم ، مشکلی پیش اومد میام سوال میپرسم باز هم، ببخشید تو زحمت افتادین هم شما و هم مهندس hadimakarem ، دست جفتتون درد نکنه ، خدا خیرتون بدهد، سر نماز و تو این شب های عزیز دعاتون میکنم.
................
و اما درباره امتحان،اتفاقا معقوله امتحانش!
آخه کلا 20 نمره رو 3 تیکه کرده تو سه تا امتحان که دو تاش رو دادیم و اخریش مونده!
اول شبیه سازی کردیم با مجموعه برنامه Aspen که 2 نمره داشت و بعدش امتحان 10 نمره ای کتبی داشتیم برای مدل سازی انواع تجهیزات مثل راکتور ها و پمپ ها و مبدل های حرارتی و.....و 26 تیر ماه هم امتحان عملی کد نویسی با متلب داریم 8 نمره ای که کد اماده با خودمون هر چند تا که خواستیم میبریم و اونجا بهمون یکی دو تا مسئله میده که اول مدل سازیش میکنیم و تعداد معلوم و مجهول و معادله ها رو بدست میاریم و بعدش چند تا از این متغیر ها و ثابت ها رو عدد هاشون رو میده و میگه حالا با متلب کد نویسیش بکنید، رانگ کوتا مرتبه 4 و چند روش دیگر مثل نصف کردن فاصله و نیوتن-رافسون و موقعیت نابجا و ..... و رو باید کد های اماده داشته باشیم که اونجا فقط در اخر جایگذاری میکنیم.
کلا برای 3 واحد 3 بار امتحان داریم میدیم !
% Runge8clear('all'), close('all'), clc
h=0.5;
x0=0;
xf = 2;
x = x0 : h : xf; x = x(:);
y0 = [4,6];
y = zeros(length(x), length(y0));
y(1,:) = y0(:).';
for i=1:length(x)-1
k1=RK(x(i),y(i,:));
k11 = k1(1);
k12 = k1(2);
ytmp1 = y(i,1) + k11*h/2;
ytmp2 = y(i,2) + k12*h/2;
ytmp = y(i,:) + k1.'*h/2;
k2=RK(x(i)+h/2,ytmp);
k21 = k2(1);
k22 = k2(2);
ytmp1 = y(i,1) + k21*h/2;
ytmp2 = y(i,2) + k22*h/2;
ytmp = y(i,:) + k2.'*h/2;
k3=RK(x(i)+h/2,ytmp);
k31 = k3(1);
k32 = k3(2);
ytmp1 = y(i,1) + k31*h/2;
ytmp2 = y(i,2) + k32*h/2;
ytmp = y(i,:) + k3.'*h;
k4=RK(x(i)+h,ytmp);
k41 = k4(1);
k42 = k4(2);
ytmp1 = y(i,1) + k41*h/2;
ytmp2 = y(i,2) + k42*h/2;
y(i+1,:)=y(i,:)+(h/6)*(k1+2*k2+2*k3+k4).';
end
xy = [x, y]
plot(x,y,'o-')
y = zeros(length(x), length(y0));
y(1,:) = y0(:).';
.'
مهندس این جوری باید بشه؟
کد:% Runge8clear('all'), close('all'), clc h=0.5; x0=0; xf = 2; x = x0 : h : xf; x = x(:); y0 = [4,6]; y = zeros(length(x), length(y0)); y(1,:) = y0(:).'; for i=1:length(x)-1 k1=RK(x(i),y(i,:)); k11 = k1(1); k12 = k1(2); ytmp1 = y(i,1) + k11*h/2; ytmp2 = y(i,2) + k12*h/2; ytmp = y(i,:) + k1.'*h/2; k2=RK(x(i)+h/2,ytmp); k21 = k2(1); k22 = k2(2); ytmp1 = y(i,1) + k21*h/2; ytmp2 = y(i,2) + k22*h/2; ytmp = y(i,:) + k2.'*h/2; k3=RK(x(i)+h/2,ytmp); k31 = k3(1); k32 = k3(2); ytmp1 = y(i,1) + k31*h/2; ytmp2 = y(i,2) + k32*h/2; ytmp = y(i,:) + k3.'*h; k4=RK(x(i)+h,ytmp); k41 = k4(1); k42 = k4(2); ytmp1 = y(i,1) + k41*h/2; ytmp2 = y(i,2) + k42*h/2; y(i+1,:)=y(i,:)+(h/6)*(k1+2*k2+2*k3+k4).'; end xy = [x, y] plot(x,y,'o-')
شرمنده اگر سوال مبتدیانه میپرسم ، آخه هنوز چند روز نمیشه دارم متلب کار میکنم و با ریزه کاری هاش آشنا نشدم!
نه دیگه؛ اون ytmp رو که به همون صورت نگه داشتید، در این صورت برنامه از ytmp1 و ytmp2 هیچ استفاده ای نمی کنه. می تونید به جاش این طور بنویسید:
کد:ytmp = [ytmp1, ytmp2];
میتونید بهم بگید این دو خط رو چطور نوشتید و چه کار میکنن برای ما این دستورات؟
خط اول یه ماتریس صفر با ابعادی که لازم داربم درست میکنه. ورودیهاش ابعاد ماتریس هستند. خط دوم آغازینه (یا شرایط اولیه) رو در سطر اول ماتریس قبلی جایگزین می کنه.کد:y = zeros(length(x), length(y0)); y(1,:) = y0(:).';
بعضی از علامت ها رو نمیدونم برای چی میگذارید اگر توضیح بدید بهم ممنون میشم از شما،مثل این علامتی که گذاشتید چند بار دقیقا چه کار میکنه ؟ چونکه تازه اولین باره دارم باهاش مواجه میشم.
این ماتریس رو ترانهاده می کنه.کد:.'
استاد به ما خیلی ساده درس داد و ظریف و منظم نبودند کد ها ولی حالا که این مدلی رو میبینم میخوام اینجوری یاد بگیرم.
هیچی اصولی کار کردن نمیشه.
در بالا با رنگ آبی جواب دادم.
% Runge8clear('all'), close('all'), clc
h=0.5;
x0=0;
xf = 2;
x = x0 : h : xf; x = x(:);
y0 = [4,6];
y = zeros(length(x), length(y0));
y(1,:) = y0(:).';
for i=1:length(x)-1
k1=RK(x(i),y(i,:));
k11 = k1(1);
k12 = k1(2);
ytmp1 = y(i,1) + k11*h/2;
ytmp2 = y(i,2) + k12*h/2;
ytmp = [ytmp1, ytmp2];
k2=RK(x(i)+h/2,ytmp);
k21 = k2(1);
k22 = k2(2);
ytmp1 = y(i,1) + k21*h/2;
ytmp2 = y(i,2) + k22*h/2;
ytmp = [ytmp1, ytmp2];
k3=RK(x(i)+h/2,ytmp);
k31 = k3(1);
k32 = k3(2);
ytmp1 = y(i,1) + k31*h/2;
ytmp2 = y(i,2) + k32*h/2;
ytmp = [ytmp1, ytmp2];
k4=RK(x(i)+h,ytmp);
k41 = k4(1);
k42 = k4(2);
ytmp1 = y(i,1) + k41*h/2;
ytmp2 = y(i,2) + k42*h/2;
ytmp = [ytmp1, ytmp2];
y(i+1,:)=y(i,:)+(h/6)*(k1+2*k2+2*k3+k4).';
end
xy = [x, y]
plot(x,y,'o-')
xy =
0 4.0000 6.0000
0.5000 3.0967 6.8647
1.0000 2.3974 7.6457
1.5000 1.8560 8.3462
2.0000 1.4368 8.9709
xy=
0 4.0000 6.0000
0.5000 3.1152 6.8577
1.0000 2.4262 7.6321
1.5000 1.8895 8.3269
2.0000 1.4716 8.9469
function E=RK4(f,x,y,h) k1=f(x,y);
k2=f(x+h/2,y+k1*h/2);
k3=f(x+h/2,y+k2*h/2);
k4=f(x+h,y+k3*h);
E=(k1+2*k2+2*k3+k4)*h/6;
clear all;close all;
clc;
h=0.5;
x0=0;
y1(1)=4;
y2(1)=6;
f=@(x,y) [-0.5*y(1);4-0.3*y(2)-0.1*y(1)];
for i=1:4
x=x0+h*i;
E=RK4(f,x,[y1(i);y2(i)],h);
y1(i+1)=y1(i)+E(1);
y2(i+1)=y2(i)+E(2);
end
H=[y1;y2]
H =
4.0000 3.1152 2.4262 1.8895 1.4716
6.0000 6.8577 7.6321 8.3269 8.9469
% Runge kutta code 4th
%dy1/dt=-0.5*y1
%dy2/dt=4-0.3*y2-0.1*y1
clear all;
close all;
clc;
% Define Function
f1=@(x,y1,y2) -0.5*y1;
f2=@(x,y1,y2) 4-0.3*y2-0.1*y1;
%Initial Conditions
x(1)=0;
y1(1)=4;
y2(1)=6;
%Step Size
h=0.5;
xf=2;
N=ceil(xf/h);
%Update Loop
for i=1:N
%Update time
x(i+1)=x(i)+h;
%Update y1 & y2
K11=f1(x(i),y1(i),y2(i));
K12=f2(x(i),y1(i),y2(i));
K21=f1(x(i)+h/2,y1(i)+h/2*K11,y2(i)+h/2*K12);
K22=f2(x(i)+h/2,y1(i)+h/2*K11,y2(i)+h/2*K12);
K31=f1(x(i)+h/2,y1(i)+h/2*K21,y2(i)+h/2*K22);
K32=f2(x(i)+h/2,y1(i)+h/2*K21,y2(i)+h/2*K22);
K41=f1(x(i)+h,y1(i)+h*K31,y2(i)+h*K32);
K42=f2(x(i)+h,y1(i)+h*K31,y2(i)+h*K32);
y1(i+1)=y1(i)+h/6*(K11+2*K21+2*K31+K41);
y2(i+1)=y2(i)+h/6*(K12+2*K22+2*K32+K42);
end
xy=[x;y1;y2]
%Plot the solution
plot(x,y1,'rp-.')
hold on
plot(x,y2,'bh:')
xlabel('x(s)')
ylabel('y1 & y2')
legend('y1','y2')
xy =
0 0.5000 1.0000 1.5000 2.0000
4.0000 3.1152 2.4262 1.8895 1.4716
6.0000 6.8577 7.6321 8.3269 8.9469
function [f1]=rk3(x,y1,y2,h);f1=-0.5*y1;
end
function [f2]=rk4(x,y1,y2)f2=4-0.3*y2-0.1*y1;
end
% Runge kutta code 4th
%dy1/dt=-0.5*y1
%dy2/dt=4-0.3*y2-0.1*y1
clear all;
close all;
clc;
%Initial Conditions
x(1)=0;
y1(1)=4;
y2(1)=6;
%Step Size
h=0.5;
xf=2;
N=ceil(xf/h);
%Update Loop
for i=1:N
%Update time
x(i+1)=x(i)+h;
%Update y1 & y2
K11=rk3(x(i),y1(i),y2(i));
K12=rk4(x(i),y1(i),y2(i));
K21=rk3(x(i)+h/2,y1(i)+h/2*K11,y2(i)+h/2*K12);
K22=rk4(x(i)+h/2,y1(i)+h/2*K11,y2(i)+h/2*K12);
K31=rk3(x(i)+h/2,y1(i)+h/2*K21,y2(i)+h/2*K22);
K32=rk4(x(i)+h/2,y1(i)+h/2*K21,y2(i)+h/2*K22);
K41=rk3(x(i)+h,y1(i)+h*K31,y2(i)+h*K32);
K42=rk4(x(i)+h,y1(i)+h*K31,y2(i)+h*K32);
y1(i+1)=y1(i)+h/6*(K11+2*K21+2*K31+K41);
y2(i+1)=y2(i)+h/6*(K12+2*K22+2*K32+K42);
end
xy=[x;y1;y2]
%Plot the solution
plot(x,y1,'rp-.')
hold on
plot(x,y2,'bh:')
xlabel('x(s)')
ylabel('y1 & y2')
legend('y1','y2')
xy =
0 0.5000 1.0000 1.5000 2.0000
4.0000 3.1152 2.4262 1.8895 1.4716
6.0000 6.8577 7.6321 8.3269 8.9469
function [f1,f2]=rk5(x,y1,y2,h);f1=-0.5*y1;
f2=4-0.3*y2-0.1*y1;
end
% Runge kutta code 4th
%dy1/dt=-0.5*y1
%dy2/dt=4-0.3*y2-0.1*y1
clear all;
close all;
clc;
%Initial Conditions
x(1)=0;
y1(1)=4;
y2(1)=6;
%Step Size
h=0.5;
xf=2;
N=ceil(xf/h);
%Update Loop
for i=1:N
%Update time
x(i+1)=x(i)+h;
%Update y1 & y2
K11=rk5(x(i),y1(i),y2(i));
K12=rk5(x(i),y1(i),y2(i));
K21=rk5(x(i)+h/2,y1(i)+h/2*K11,y2(i)+h/2*K12);
K22=rk5(x(i)+h/2,y1(i)+h/2*K11,y2(i)+h/2*K12);
K31=rk5(x(i)+h/2,y1(i)+h/2*K21,y2(i)+h/2*K22);
K32=rk5(x(i)+h/2,y1(i)+h/2*K21,y2(i)+h/2*K22);
K41=rk5(x(i)+h,y1(i)+h*K31,y2(i)+h*K32);
K42=rk5(x(i)+h,y1(i)+h*K31,y2(i)+h*K32);
y1(i+1)=y1(i)+h/6*(K11+2*K21+2*K31+K41);
y2(i+1)=y2(i)+h/6*(K12+2*K22+2*K32+K42);
end
xy=[x;y1;y2]
%Plot the solution
plot(x,y1,'rp-.')
hold on
plot(x,y2,'bh:')
xlabel('x(s)')
ylabel('y1 & y2')
legend('y1','y2')
xy =
0 0.5000 1.0000 1.5000 2.0000
4.0000 3.1152 2.4262 1.8895 1.4716
6.0000 5.1152 4.4262 3.8895 3.4716
سلام مهندس جان.
ممنون از راهنمایی.
بنده این چند روز داشتم چند تا فایل آموزشی میخوندم.
اون کاری رو که گفتید با کدی که نوشتین انجام دادم ولی جوابش فرق کرد!!!
کد:% Runge8clear('all'), close('all'), clc h=0.5; x0=0; xf = 2; x = x0 : h : xf; x = x(:); y0 = [4,6]; y = zeros(length(x), length(y0)); y(1,:) = y0(:).'; for i=1:length(x)-1 k1=RK(x(i),y(i,:)); k11 = k1(1); k12 = k1(2); ytmp1 = y(i,1) + k11*h/2; ytmp2 = y(i,2) + k12*h/2; ytmp = [ytmp1, ytmp2]; k2=RK(x(i)+h/2,ytmp); k21 = k2(1); k22 = k2(2); ytmp1 = y(i,1) + k21*h/2; ytmp2 = y(i,2) + k22*h/2; ytmp = [ytmp1, ytmp2]; k3=RK(x(i)+h/2,ytmp); k31 = k3(1); k32 = k3(2); ytmp1 = y(i,1) + k31*h/2; ytmp2 = y(i,2) + k32*h/2; ytmp = [ytmp1, ytmp2]; k4=RK(x(i)+h,ytmp); k41 = k4(1); k42 = k4(2); ytmp1 = y(i,1) + k41*h/2; ytmp2 = y(i,2) + k42*h/2; ytmp = [ytmp1, ytmp2]; y(i+1,:)=y(i,:)+(h/6)*(k1+2*k2+2*k3+k4).'; end xy = [x, y] plot(x,y,'o-')
این جوابی هست که بدست اوردم :
کد:xy = 0 4.0000 6.0000 0.5000 3.0967 6.8647 1.0000 2.3974 7.6457 1.5000 1.8560 8.3462 2.0000 1.4368 8.9709
ولی کد اصلی جواب درست رو میده که این هستش:
کد:xy= 0 4.0000 6.0000 0.5000 3.1152 6.8577 1.0000 2.4262 7.6321 1.5000 1.8895 8.3269 2.0000 1.4716 8.9469
خب ببینید مهندس جان مشکل من اینه که من الان این مثال رو به چند روش حل کردم ولی این با همه شون فرق داره!
مدل کلیش نه! مدل نوشتنش!
کاملا حرفه ای هستش کد شما ولی من نمیتونم این رو بفهمم که چرا شما این مدلی نوشتین که ماتریس رو ترانهاده میکنید!
میتونیدبهم توضیح بدین که دقیقا چرا این کار رو کردین؟
.................
اینجا میخوام برای شما چند مدل دیگه ای که نوشتم رو بزارم ،البته مبتدیانه نوشته شدن :
این کد اولی همون راه آسونه هستش :
کد:function E=RK4(f,x,y,h) k1=f(x,y); k2=f(x+h/2,y+k1*h/2); k3=f(x+h/2,y+k2*h/2); k4=f(x+h,y+k3*h); E=(k1+2*k2+2*k3+k4)*h/6;
و
کد:clear all;close all; clc; h=0.5; x0=0; y1(1)=4; y2(1)=6; f=@(x,y) [-0.5*y(1);4-0.3*y(2)-0.1*y(1)]; for i=1:4 x=x0+h*i; E=RK4(f,x,[y1(i);y2(i)],h); y1(i+1)=y1(i)+E(1); y2(i+1)=y2(i)+E(2); end H=[y1;y2]
جوابش هم درسته :
کد:H = 4.0000 3.1152 2.4262 1.8895 1.4716 6.0000 6.8577 7.6321 8.3269 8.9469
ولی خب من این مدلیش رو نمیخوام!
(یه مشکلی که داره اینکه نمیتونم x ها رو وارد ماتریس بکنم! ارور میده!)
.........
مدل بعدی این هستش که همه چی رو تو یه اسکریپت نوشتم!
کد:% Runge kutta code 4th %dy1/dt=-0.5*y1 %dy2/dt=4-0.3*y2-0.1*y1 clear all; close all; clc; % Define Function f1=@(x,y1,y2) -0.5*y1; f2=@(x,y1,y2) 4-0.3*y2-0.1*y1; %Initial Conditions x(1)=0; y1(1)=4; y2(1)=6; %Step Size h=0.5; xf=2; N=ceil(xf/h); %Update Loop for i=1:N %Update time x(i+1)=x(i)+h; %Update y1 & y2 K11=f1(x(i),y1(i),y2(i)); K12=f2(x(i),y1(i),y2(i)); K21=f1(x(i)+h/2,y1(i)+h/2*K11,y2(i)+h/2*K12); K22=f2(x(i)+h/2,y1(i)+h/2*K11,y2(i)+h/2*K12); K31=f1(x(i)+h/2,y1(i)+h/2*K21,y2(i)+h/2*K22); K32=f2(x(i)+h/2,y1(i)+h/2*K21,y2(i)+h/2*K22); K41=f1(x(i)+h,y1(i)+h*K31,y2(i)+h*K32); K42=f2(x(i)+h,y1(i)+h*K31,y2(i)+h*K32); y1(i+1)=y1(i)+h/6*(K11+2*K21+2*K31+K41); y2(i+1)=y2(i)+h/6*(K12+2*K22+2*K32+K42); end xy=[x;y1;y2] %Plot the solution plot(x,y1,'rp-.') hold on plot(x,y2,'bh:') xlabel('x(s)') ylabel('y1 & y2') legend('y1','y2')
جوابش هم درسته :
کد:xy = 0 0.5000 1.0000 1.5000 2.0000 4.0000 3.1152 2.4262 1.8895 1.4716 6.0000 6.8577 7.6321 8.3269 8.9469
ولی خب نتونستم یه function درست کنم برای معادلات! چونکه همه اش ارور های مسخره میده!
.......................
یه مدل دیگه نوشتم که function ها رو دو قسمت کردم و جواب داد.
این کدش:
وکد:function [f1]=rk3(x,y1,y2,h);f1=-0.5*y1; end
وکد:function [f2]=rk4(x,y1,y2)f2=4-0.3*y2-0.1*y1; end
کد:% Runge kutta code 4th %dy1/dt=-0.5*y1 %dy2/dt=4-0.3*y2-0.1*y1 clear all; close all; clc; %Initial Conditions x(1)=0; y1(1)=4; y2(1)=6; %Step Size h=0.5; xf=2; N=ceil(xf/h); %Update Loop for i=1:N %Update time x(i+1)=x(i)+h; %Update y1 & y2 K11=rk3(x(i),y1(i),y2(i)); K12=rk4(x(i),y1(i),y2(i)); K21=rk3(x(i)+h/2,y1(i)+h/2*K11,y2(i)+h/2*K12); K22=rk4(x(i)+h/2,y1(i)+h/2*K11,y2(i)+h/2*K12); K31=rk3(x(i)+h/2,y1(i)+h/2*K21,y2(i)+h/2*K22); K32=rk4(x(i)+h/2,y1(i)+h/2*K21,y2(i)+h/2*K22); K41=rk3(x(i)+h,y1(i)+h*K31,y2(i)+h*K32); K42=rk4(x(i)+h,y1(i)+h*K31,y2(i)+h*K32); y1(i+1)=y1(i)+h/6*(K11+2*K21+2*K31+K41); y2(i+1)=y2(i)+h/6*(K12+2*K22+2*K32+K42); end xy=[x;y1;y2] %Plot the solution plot(x,y1,'rp-.') hold on plot(x,y2,'bh:') xlabel('x(s)') ylabel('y1 & y2') legend('y1','y2')
جوابش هم درست در میاد:
کد:xy = 0 0.5000 1.0000 1.5000 2.0000 4.0000 3.1152 2.4262 1.8895 1.4716 6.0000 6.8577 7.6321 8.3269 8.9469
ولی باز این مدلی هم نمیخوام.
..............................
این هم یه مدل دیگه که من میخوام کدم به این شکل باشه با جواب درست!
میخوام به جای 2 تا function یکی بنویسم .
این هم کدش :
وکد:function [f1,f2]=rk5(x,y1,y2,h);f1=-0.5*y1; f2=4-0.3*y2-0.1*y1; end
کد:% Runge kutta code 4th %dy1/dt=-0.5*y1 %dy2/dt=4-0.3*y2-0.1*y1 clear all; close all; clc; %Initial Conditions x(1)=0; y1(1)=4; y2(1)=6; %Step Size h=0.5; xf=2; N=ceil(xf/h); %Update Loop for i=1:N %Update time x(i+1)=x(i)+h; %Update y1 & y2 K11=rk5(x(i),y1(i),y2(i)); K12=rk5(x(i),y1(i),y2(i)); K21=rk5(x(i)+h/2,y1(i)+h/2*K11,y2(i)+h/2*K12); K22=rk5(x(i)+h/2,y1(i)+h/2*K11,y2(i)+h/2*K12); K31=rk5(x(i)+h/2,y1(i)+h/2*K21,y2(i)+h/2*K22); K32=rk5(x(i)+h/2,y1(i)+h/2*K21,y2(i)+h/2*K22); K41=rk5(x(i)+h,y1(i)+h*K31,y2(i)+h*K32); K42=rk5(x(i)+h,y1(i)+h*K31,y2(i)+h*K32); y1(i+1)=y1(i)+h/6*(K11+2*K21+2*K31+K41); y2(i+1)=y2(i)+h/6*(K12+2*K22+2*K32+K42); end xy=[x;y1;y2] %Plot the solution plot(x,y1,'rp-.') hold on plot(x,y2,'bh:') xlabel('x(s)') ylabel('y1 & y2') legend('y1','y2')
این هم جوابش :
کد:xy = 0 0.5000 1.0000 1.5000 2.0000 4.0000 3.1152 2.4262 1.8895 1.4716 6.0000 5.1152 4.4262 3.8895 3.4716
ولی جوابش غلطه کلا برای y2 به جای زیاد شدن کم میشه عددش!
مثل اینکه به درستی از معادله دومی استتفاده نمیکنه با اینکه معادله اولی خرابش میکنه! نمیدونم والله!شما اگر بتونید مشکل رو بهم بگید ممنونتون میشم.
شرمنده همه عزیزان و مخصوصا مهندس عزیز اگر پست بسیار طولانی شد.
% Runge8sc
clear('all'), close('all'), clc
h=0.5;
x0=0;
xf = 2;
x = x0 : h : xf; x = x(:);
y0 = [4,6];
y = zeros(length(x), length(y0));
y(1,:) = y0(:).';
for i=1:length(x)-1
k1=RK(x(i),y(i,:));
k11 = k1(1);
k12 = k1(2);
ytmp1 = y(i,1) + k11*h/2;
ytmp2 = y(i,2) + k12*h/2;
ytmp = [ytmp1, ytmp2];
k2=RK(x(i)+h/2,ytmp);
k21 = k2(1);
k22 = k2(2);
ytmp1 = y(i,1) + k21*h/2;
ytmp2 = y(i,2) + k22*h/2;
ytmp = [ytmp1, ytmp2];
k3=RK(x(i)+h/2,ytmp);
k31 = k3(1);
k32 = k3(2);
ytmp1 = y(i,1) + k31*h;
ytmp2 = y(i,2) + k32*h;
ytmp = [ytmp1, ytmp2];
k4=RK(x(i)+h,ytmp);
k41 = k4(1);
k42 = k4(2);
ytmp1 = y(i,1) + (h/6)*(k11+2*k21+2*k31+k41);
ytmp2 = y(i,2) + (h/6)*(k12+2*k22+2*k32+k42);
y(i+1,:) = [ytmp1, ytmp2];
end
xy = [x, y]
plot(x,y,'o-')
سلام،
من فقط اون دومی رو برات اصلاح کردم؛ دو جا کم دقتی کردید. باقیش من الآن وقت ندارم چک کنم.
کد:% Runge8sc clear('all'), close('all'), clc h=0.5; x0=0; xf = 2; x = x0 : h : xf; x = x(:); y0 = [4,6]; y = zeros(length(x), length(y0)); y(1,:) = y0(:).'; for i=1:length(x)-1 k1=RK(x(i),y(i,:)); k11 = k1(1); k12 = k1(2); ytmp1 = y(i,1) + k11*h/2; ytmp2 = y(i,2) + k12*h/2; ytmp = [ytmp1, ytmp2]; k2=RK(x(i)+h/2,ytmp); k21 = k2(1); k22 = k2(2); ytmp1 = y(i,1) + k21*h/2; ytmp2 = y(i,2) + k22*h/2; ytmp = [ytmp1, ytmp2]; k3=RK(x(i)+h/2,ytmp); k31 = k3(1); k32 = k3(2); ytmp1 = y(i,1) + k31*h; ytmp2 = y(i,2) + k32*h; ytmp = [ytmp1, ytmp2]; k4=RK(x(i)+h,ytmp); k41 = k4(1); k42 = k4(2); ytmp1 = y(i,1) + (h/6)*(k11+2*k21+2*k31+k41); ytmp2 = y(i,2) + (h/6)*(k12+2*k22+2*k32+k42); y(i+1,:) = [ytmp1, ytmp2]; end xy = [x, y] plot(x,y,'o-')
function f = RK8(t,CA)TAU=2;
K=0.5;
CA0=1.8;
f = [1/TAU*(CA0-CA(1))-K*CA(1);1/TAU*(CA(1)-CA(2))-K*CA(2);1/TAU*(CA(2)-CA(3))-K*CA(3)];
clear('all'), close('all'), clc
h=0.1;
t0=0;
tf = 2.9;
t = t0 : h : tf; t = t(:);
CA1 = [0.4,0.2,0.1];
CA = zeros(length(t), length(CA1));
CA(1,:) = CA1(:).';
N=ceil(tf/h);
for i=1:N
k1=RK8(t(i),CA(i,:));
k11 = k1(1);
k12 = k1(2);
k13 = k1(3);
CAtmp1 = CA(i,1) + k11*h/2;
CAtmp2 = CA(i,2) + k12*h/2;
CAtmp3 = CA(i,3) + k13*h/2;
CAtmp = [CAtmp1, CAtmp2,CAtmp3];
k2=RK8(t(i)+h/2,CAtmp);
k21 = k2(1);
k22 = k2(2);
k23 = k2(3);
CAtmp1 = CA(i,1) + k21*h/2;
CAtmp2 = CA(i,2) + k22*h/2;
CAtmp3 = CA(i,3) + k23*h/2;
CAtmp = [CAtmp1, CAtmp2,CAtmp3];
k3=RK8(t(i)+h/2,CAtmp);
k31 = k3(1);
k32 = k3(2);
k33 = k3(3);
CAtmp1 = CA(i,1) + k31*h;
CAtmp2 = CA(i,2) + k32*h;
CAtmp3 = CA(i,3) + k33*h;
CAtmp = [CAtmp1, CAtmp2,CAtmp3];
k4=RK8(t(i)+h,CAtmp);
k41 = k4(1);
k42 = k4(2);
k43 = k4(3);
CAtmp1 = CA(i,1) + (h/6)*(k11+2*k21+2*k31+k41);
CAtmp2 = CA(i,2) + (h/6)*(k12+2*k22+2*k32+k42);
CAtmp3 = CA(i,3) + (h/6)*(k13+2*k23+2*k33+k43);
CA(i+1,:) = [CAtmp1, CAtmp2,CAtmp3];
end
tCA = [t, CA]
plot(t,CA,'o-')
tCA =
0 0.400000000000000 0.200000000000000 0.100000000000000
0.100000000000000 0.447581250000000 0.201169791666667 0.100019270833333
0.200000000000000 0.490634549296875 0.204380918085938 0.100143457419705
0.300000000000000 0.529590788999411 0.209234268703024 0.100449806855966
0.400000000000000 0.564839855541255 0.215388207553221 0.100990646230996
0.500000000000000 0.596734532788310 0.222551248716383 0.101798307935674
0.600000000000000 0.625594032811842 0.230475605331308 0.102889259307639
0.700000000000000 0.651707190664386 0.238951512240097 0.104267551429592
0.800000000000000 0.675335355132786 0.247802233440856 0.105927686492374
0.900000000000000 0.696715004399962 0.256879675426173 0.107856989691625
1.000000000000000 0.716060112793751 0.266060536304158 0.110037559907895
1.100000000000000 0.733564292310016 0.275242928454424 0.112447863212385
1.200000000000000 0.749402730343064 0.284343419466639 0.115064024356605
1.300000000000000 0.763733943016792 0.293294442336916 0.117860863679840
1.400000000000000 0.776701361664456 0.302042031440014 0.120812720158337
1.500000000000000 0.788434768335063 0.310543845727073 0.123894095497385
1.600000000000000 0.799051594693377 0.318767444985496 0.127080149120063
1.700000000000000 0.808658097313369 0.326688788898650 0.130347069536140
1.800000000000000 0.817350421127785 0.334290932111034 0.133672343795083
1.900000000000000 0.825215561677212 0.341562891586541 0.137034943462841
2.000000000000000 0.832332235789105 0.348498665285483 0.140415442746567
2.100000000000000 0.838771669400824 0.355096383617746 0.143796081966557
2.200000000000000 0.844598310411468 0.361357577288197 0.147160787489597
2.300000000000000 0.849870473696937 0.367286547066375 0.150495157447918
2.400000000000000 0.854640924743752 0.372889822712010 0.153786421034178
2.500000000000000 0.858957407742825 0.378175699795127 0.157023377851159
2.600000000000000 0.862863123428498 0.383153844485338 0.160196322675471
2.700000000000000 0.866397161445234 0.387834957568613 0.163296960041503
2.800000000000000 0.869594891569202 0.392230489998041 0.166318312242438
2.900000000000000 0.872488317700248 0.396352403213013 0.169254623659750
این مثال بالا رو با دو روش دیگه حل کردم و همین جواب رو داد پس جوابش درسته فکر کنم.
............................
الان دارم روی دو تا مثال دیگه کار میکنم ولی روش این مثال ها رانگ کوتا نیستش.
این دو تا مثال از طریق روش نصف کردن و روش نیوتن-رافسون بدست میان که من تا به حال مدلی مشابه شون حل نکردم که بدونم چجوری حل میشن و این زبان فورتن هم رو هم نمیفهمم که چی میگه!
میتونید این دو تا کدی که با فورتن نوشته شدن رو به صورت کد برای متلب بنویسید؟
بنده کل اون صفحاتی که مورد نیاز این دو مثال هست اعم از متن توضیح دو روش و صورت سوال و کد و جواب هاشون رو عکس گرفتم از کتاب و تو یه پوشه قرار دادم.
ممنونتون میشم.
http://s7.picofile.com/file/8259049450/2Question.rar.html
clear all;close all;
clc;
U = 100;
A = 40;
TC1 = 35;
TH1 = 300;
TC2 = 125;
Cp_warm = 1000;
Cp_cold = 4197;
m_warm = 2;
Error = 10;
TH2 = 150;
while Error>0.001
Q2 = Cp_warm*m_warm*(TH1-TH2) ;
m_cold = (Q2/(Cp_cold*(TC2-TC1)));
T_LM = ((TH1-TC2)-(TH2-TC1))/(log((TH1-TC2)/(TH2-TC1)));
Q1 = U*A*T_LM;
TH2_new = TH1-(Q1/(Cp_warm*m_warm));
Error = abs(TH2_new - TH2);
TH2 = TH2_new;
end
disp(TH2);
disp(m_cold);
Thread starter | عنوان | تالار | پاسخ ها | تاریخ |
---|---|---|---|---|
F | سوالات و مشکلات در متلب | MATLAB | 1 | |
N | راهنمایی پیاده سازی کلاسبند جنگل تصادفی و خوشه بند DBSCAN در متلب | MATLAB | 0 | |
ب | کمک در برنامه نویسی متلب...فوری | MATLAB | 0 | |
A | شبیه سازی متغییر گسسته وپیوسته باهم برای یک مسئله در متلب 2014 | MATLAB | 0 | |
A | درخواست کمک در پروژه متلب.فوری | MATLAB | 4 |