(هادی)
کاربر فعال تالار ریاضی ,
احتمالا خطای عددی باشه. با عوض کردن دستور حل معادله دیفرانسیل، یا با کوچک کردن گام زمانی حل، باید درست بشه
سئوالات و مشکلات در متلب (MATLAB)
- شروع کننده موضوع arashonline
- تاریخ شروع
از روش ضمنی استفاده کردم و آخرt,c) plot) را نوشتم.از ode استفاده نکردم
توی اون پست قبلی هم گفتم از روش صریح استفاده کنید.
اگه می خواید حتماً از روش ضمنی استفاده کنید، چند تا کار که در CFD مرسومه می تونید انجام بدید:
1. با استفاده از تحلیل پایداری (مثلاً روش نیومن) بازه یا مقدار مناسب گامها رو پیدا کنید. پیشنهاد من اینه که به جای این کار برنامه رو با جایگشتهای مختلف گامهای dr، dz و dt اجرا کنید و هر جوابی که بهتر بود رو نگه دارید. توجه کنید که در اینگونه مسائل کوچک کردن یک گام لزوماً منجر به جواب بهتری نمی شه؛ ممکنه بزرگتر کردن اون گام منجر به جواب بهتری بشه.
2. اگه همه این کارها رو انجام دادید و نتیجه مطلوب رو به دست نیاوردید، می تونید جواب رو از یک فیلتر low pass عبور بدید؛ البته اگه این کار در حین حل معادلات انجام بشه بهتره تا بعد از حل معادلات.
dt .dz.dr را تغییر دادم.هرچی کمتر باشن جوابا بهترن...ببخشید low pass چی هست؟چجوری بکار میره؟
با سلام .
نماز و روزه و طاعات و عبادات شما عزیزان قبول درگاه حق انشاالله.
میخواستم این کد رو برام تصحیح بکنید که بتونم اعداد تو جدول رو بدست بیارم،
یه دستگاه دو معادله ای هستش که باید با روش رانگ کوتا مرتبه 4 حل بشه که دو مدل داره یکی اینکه رانگ کوتا رو به صورت یه تابع تعریف بکنیم و نو برنامه ای که مینویسیم فقط معادلات رو جایگذرای میکنیم و یکی دیگه اینکه تو برنامه اصلی معادلات رو مینویسیم و تو برنامه ای که داریم تمام k ها رو حساب میکنیم.
من الان این مدل دومی رو نیاز دارم.
نمیخوام پیشرفته باشه، مثل همین که خودم نوشتم فقط تصحیح بشه .
با تشکر فراوان.
http://s7.picofile.com/file/8257314892/kode_2_moadeleie.rar.html
نماز و روزه و طاعات و عبادات شما عزیزان قبول درگاه حق انشاالله.
میخواستم این کد رو برام تصحیح بکنید که بتونم اعداد تو جدول رو بدست بیارم،
یه دستگاه دو معادله ای هستش که باید با روش رانگ کوتا مرتبه 4 حل بشه که دو مدل داره یکی اینکه رانگ کوتا رو به صورت یه تابع تعریف بکنیم و نو برنامه ای که مینویسیم فقط معادلات رو جایگذرای میکنیم و یکی دیگه اینکه تو برنامه اصلی معادلات رو مینویسیم و تو برنامه ای که داریم تمام k ها رو حساب میکنیم.
من الان این مدل دومی رو نیاز دارم.
نمیخوام پیشرفته باشه، مثل همین که خودم نوشتم فقط تصحیح بشه .
با تشکر فراوان.
http://s7.picofile.com/file/8257314892/kode_2_moadeleie.rar.html
آخرین ویرایش:
am_salmani
عضو جدید
سلام ممکنه این لینک به دردتون بخوره
آخرین ویرایش توسط مدیر:
(هادی)
کاربر فعال تالار ریاضی ,
سلام
نماز و روزه شما هم قبول باشه
به یک همچین فرمی باید بنویسید:
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)];
موفق باشید
خیلی ممنون مهندس ، حل شد، ولی باز اون چیزی که میخوام نشد!
من نیاز دارم که 8 تا k داشته باشم و y ها به صورت 1 و 2 باشه.
چونکه قراره در اخر 40 تا k بنویسم برای 10 معادله ای.
اون روشی که تو function خود کد رانگ کوتا رو مینویسیم خیلی بهتره چونکه اون هم فقط تعریف اصلی رانگ کوتا رو داریم و فقط تو برنامه ای که مینویسیم برای هر معادله یه خط مینویسیم ولی استاد میخواد که ما این مدلی که سخت تر هستش رو یاد بگیرم که برای هر معادله 4 تا k و جواب آخر بنویسیم.
این روش ساده هستش که خیلی هم مختصر و مفید و تمیز هستش ولی استاد قبول نمیکنه :
.................
function E=RK4(f,x,y,h)
k1=f(x,y);
k2=f(x+h/2,y+k1*h/2);
k3=f(x+h/4,y+k2*h/2);
k4=f(x+h/2,y+k3*h);
E=h*(k1+2*k2+2*k3+k4)/6;
.....................
clc;
clear all;
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
disp([y1;y2])
..............
آخرین ویرایش:
مهندس یعنی باید 2 تا function درست بکنم و داخل هر کدوم جداگونه تابع ها رو تعریف بکنم؟
باز برمیگردم به خونه اول فکر کنم!
من همون روز که بهمون گفت چی میخواد از ما 2 تا function نوشتم و یه فایل اسکریپت و وقتی میخواستم 8 تا k ها رو حساب بکنم برای هر 4 تاش به یکیشون لینک میدادم ولی استاد گفت نه این خیلی مبتدیانه هستش و همه اش باید تو یک function باشه!
................
یا اینکه منظورتون اینه که این 2 تا f1 و f2 رو تو خود فایل اسکریپتی که دارم مینویسم تعریف بکنم و نه داخل function ؟
(هادی)
کاربر فعال تالار ریاضی ,
بالاخره یا باید متغیرها رو از هم جدا کنی و هر کدوم رو جداگانه پهاربار صدا بزنی، یا همه رو یک کاسه کنی توی یک تابع
اگه همه رو یکی کردی (که منطقی تر هم همینه)، همون تابع رو چهار بار صدا میزنی ...
به عقلم جور در نمیاد که همه رو بریزی توی یک تابع، ولی به جای چهار بار، هشت بار صدا بزنی اون تابع رو
والله به خدا به عقل من هم جور در نمیاد !
ببین الان من اینکار رو کردم که 2 تا function جدا گونه ساختم و تو کدی که مینویسم برای هر 4 تا k یکی از function ها رو انتخاب میکنیم
....................
function [f1]=RK6(x,y,t,h);
f1 = [-0.5*y];
....................
function [f2]=RK7(x,y,t,h);
f2 = [4-0.3*t-0.1*y];
...................
clear all
clc
h=0.5;
x0=0;
y=4;
t=6;
for i=1:4;
x=x0+i*h;
k1=RK6(x,y,t);
k2=RK6(x+0.5*h,y+0.5*h*k1,t+0.5*h*k1);
k3=RK6(x+0.5*h,y+0.5*h*k2,t+0.5*h*k2);
k4=RK6(x+h,y+k3*h,t+k3*h);
k5=RK7(x,y,t);
k6=RK7(x+0.5*h,y+0.5*h*k5,t+0.5*h*k5);
k7=RK7(x+0.5*h,y+0.5*h*k6,t+0.5*h*k6);
k8=RK7(x+h,y+k7*h,t+k7*h);
y=y+(1/6)*(k1+2*k2+2*k3+k4)*h
t=t+(1/6)*(k5+2*k6+2*k7+k8)*h
end
...................
جواب بدست اومد بدون ارور دادن ولی جواب مثل قبلی ها درست در نمیاد! یه اختلاف کوچکی داره که نمیدونم به خاطر چی هستش!
y =
3.1152
t =
6.8157
y =
2.4262
t =
7.5606
y =
1.8895
t =
8.2355
y =
1.4716
t =
8.8429
......................
این رو باید بدست بیارم که با اون راه حل بدست اومده :
4.0000 3.1152 2.4262 1.8895 1.4716
6.0000 6.8577 7.6321 8.3269 8.9469
.................
به جای y1,y2, گذاشتم y,t میخواستم سریع بنویسم.
آخرین ویرایش:
کد:
>> help sgolayانجام دادم ولی جواب اصلا تغییری نکرد!
درست میگید برای f1 فقط معادله اش y داره که الان t هاش رو حذف کردم که کلا جوابش تغییری نکرد!
کلا y ها یا همون y1 درست در میاد!
مشکل t یا همون y2 هستش ! برای f2 چونکه هم y1,y2 یا y,t داره نمیتونم y رو حذف بکنم!
آخرین ویرایش:
در مورد روزه؛ از 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 و 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 = [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 به جای زیاد شدن کم میشه عددش!
مثل اینکه به درستی از معادله دومی استتفاده نمیکنه با اینکه معادله اولی خرابش میکنه! نمیدونم والله!شما اگر بتونید مشکل رو بهم بگید ممنونتون میشم.
شرمنده همه عزیزان و مخصوصا مهندس عزیز اگر پست بسیار طولانی شد.
آخرین ویرایش:
سلام مهندس جان.
ممنون از راهنمایی.
بنده این چند روز داشتم چند تا فایل آموزشی میخوندم.
اون کاری رو که گفتید با کدی که نوشتین انجام دادم ولی جوابش فرق کرد!!!
کد:% 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-')
آخرین ویرایش:
محاسبات انتهای حلقه رو هم به شکل اسکالر درآوردم که مشابه بقیه بخشها بشه.
سلام مهندس جان.
عید شما مبارک باشه انشاالله.
خیلی ممنون به خاطر اصلاح کد و این اصلاح دومی خیلی بهترش کرد که انتهای حلقه رو اسکالرش کردین.
اگر وقت کردین آخرین کد پست قبلیم رو هم یه نگاه بکنید خیلی عالی میشه و سپاسگذار شما میشم.
.........
مهندس این کد رو میتونید چک بکنید ببینید درسته یا خیر؟
این صورت سوالشه :
http://s6.uplod.ir/i/00797/yu3dtcxlefsb.jpg
این عکس ها هم مربوط میشن به کتاب اصلی و نوشتن کد به صورت فورتن( درسته؟) و نتایج حاصل از اون:
http://s6.uplod.ir/i/00797/0ee44dbhovbu.png
http://s6.uplod.ir/i/00797/r3mcpx8zpoml.png
این function :
این هم Script :
این هم جوابش :
این هم plot :
http://s6.uplod.ir/i/00797/tp4q77lhdhi2.jpg
عید شما مبارک باشه انشاالله.
خیلی ممنون به خاطر اصلاح کد و این اصلاح دومی خیلی بهترش کرد که انتهای حلقه رو اسکالرش کردین.
اگر وقت کردین آخرین کد پست قبلیم رو هم یه نگاه بکنید خیلی عالی میشه و سپاسگذار شما میشم.
.........
مهندس این کد رو میتونید چک بکنید ببینید درسته یا خیر؟
این صورت سوالشه :
http://s6.uplod.ir/i/00797/yu3dtcxlefsb.jpg
این عکس ها هم مربوط میشن به کتاب اصلی و نوشتن کد به صورت فورتن( درسته؟) و نتایج حاصل از اون:
http://s6.uplod.ir/i/00797/0ee44dbhovbu.png
http://s6.uplod.ir/i/00797/r3mcpx8zpoml.png
این function :
کد:
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)];این هم Script :
کد:
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این هم plot :
http://s6.uplod.ir/i/00797/tp4q77lhdhi2.jpg
آخرین ویرایش:
این مثال بالا رو با دو روش دیگه حل کردم و همین جواب رو داد پس جوابش درسته فکر کنم.
............................
الان دارم روی دو تا مثال دیگه کار میکنم ولی روش این مثال ها رانگ کوتا نیستش.
این دو تا مثال از طریق روش نصف کردن و روش نیوتن-رافسون بدست میان که من تا به حال مدلی مشابه شون حل نکردم که بدونم چجوری حل میشن و این زبان فورتن هم رو هم نمیفهمم که چی میگه!
میتونید این دو تا کدی که با فورتن نوشته شدن رو به صورت کد برای متلب بنویسید؟
بنده کل اون صفحاتی که مورد نیاز این دو مثال هست اعم از متن توضیح دو روش و صورت سوال و کد و جواب هاشون رو عکس گرفتم از کتاب و تو یه پوشه قرار دادم.
ممنونتون میشم.
http://s7.picofile.com/file/8259049450/2Question.rar.html
............................
الان دارم روی دو تا مثال دیگه کار میکنم ولی روش این مثال ها رانگ کوتا نیستش.
این دو تا مثال از طریق روش نصف کردن و روش نیوتن-رافسون بدست میان که من تا به حال مدلی مشابه شون حل نکردم که بدونم چجوری حل میشن و این زبان فورتن هم رو هم نمیفهمم که چی میگه!
میتونید این دو تا کدی که با فورتن نوشته شدن رو به صورت کد برای متلب بنویسید؟
بنده کل اون صفحاتی که مورد نیاز این دو مثال هست اعم از متن توضیح دو روش و صورت سوال و کد و جواب هاشون رو عکس گرفتم از کتاب و تو یه پوشه قرار دادم.
ممنونتون میشم.
http://s7.picofile.com/file/8259049450/2Question.rar.html
سلام،
ظاهراً پیغامی که برام گذاشته بودید در مورد این پست هست. من الآن سرم خیلی شلوغه، اما این روشهایی که اسم بردید، دو تا روش حل عددی معادله های جبری (اسکالر) هستش که توی همه کتابهای محاسبات عددی هست. توی کتاب من هم هست. برنامه هاش رو می تونید از لینک امضام دانلود کنید. بعد از unzip کردن برنامه ها، اسم تابعهایی که شما می خواید RootBis و FNewton هست.
سلام به دوستان عزیز میخواستم بدونم میتونید این کد رو برام تصحیح بکنید؟
مربوط میشه به یه مبدل حرارتی غیر هم جهت پوسته لوله ای که با روش LMTD حل میشه!
مشخصاتش اینه :
دمای ورودی سیال داغ: ۳۰۰
دبی جرمی: ۲
ظرفیت گرمایی: ۱۰۰۰
دمای ورودی سیال سرد: ۳۵
دمای خروجی سیال سرد: ۱۲۵
ظرفیت گرمایی: ۴۱۹۷
مساحت مبدل: ۴۰
ضریب انتقال حرارت کلی: ۱۰۰
چیزی که میخواییم دبی جرمی سیال سرد هستش!
دما ها به سانتی گراد هستن!
جواب من واگرا بدست میاد!
همین جور Busy میره تا چند ساعت!
باید اول دمای دوم سیال گرم رو حدس بزنیم و بعدش دبی جرمی سیال سرد رو بدست بیاریم.
این کد بنده هستش :
مربوط میشه به یه مبدل حرارتی غیر هم جهت پوسته لوله ای که با روش LMTD حل میشه!
مشخصاتش اینه :
دمای ورودی سیال داغ: ۳۰۰
دبی جرمی: ۲
ظرفیت گرمایی: ۱۰۰۰
دمای ورودی سیال سرد: ۳۵
دمای خروجی سیال سرد: ۱۲۵
ظرفیت گرمایی: ۴۱۹۷
مساحت مبدل: ۴۰
ضریب انتقال حرارت کلی: ۱۰۰
چیزی که میخواییم دبی جرمی سیال سرد هستش!
دما ها به سانتی گراد هستن!
جواب من واگرا بدست میاد!
همین جور Busy میره تا چند ساعت!
باید اول دمای دوم سیال گرم رو حدس بزنیم و بعدش دبی جرمی سیال سرد رو بدست بیاریم.
این کد بنده هستش :
کد:
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);
آخرین ویرایش:
Similar threads
| Thread starter | عنوان | تالار | پاسخ ها | تاریخ |
|---|---|---|---|---|
| F | سوالات و مشکلات در متلب | MATLAB | 1 | |
| N | راهنمایی پیاده سازی کلاسبند جنگل تصادفی و خوشه بند DBSCAN در متلب | MATLAB | 0 | |
| ب | کمک در برنامه نویسی متلب...فوری | MATLAB | 0 | |
| A | شبیه سازی متغییر گسسته وپیوسته باهم برای یک مسئله در متلب 2014 | MATLAB | 0 | |
| A | درخواست کمک در پروژه متلب.فوری | MATLAB | 4 |
Similar threads
-
-
راهنمایی پیاده سازی کلاسبند جنگل تصادفی و خوشه بند DBSCAN در متلب
- شروع شده توسط nana68
- پاسخ ها: 0
-
-
شبیه سازی متغییر گسسته وپیوسته باهم برای یک مسئله در متلب 2014
- شروع شده توسط a1b2c3650917
- پاسخ ها: 0
-