السلام عليكم ورحمة الله وبركاته… طلبتكم يا اخوة وبليييز لا تردوني محد يعلم بحالي غير رب العالمين
عندي بروجكت وتقريبا مخلص اكثر من النص بس باقيلي شغلات بسيطة معكرة مزاجي والله مو مخليتني انام الليل…
الماتلاب كود اللي سويته يشتغل ويطلع لي graphs وكل شي تمام وحلو… بس المشكلة اني اواجه بعض errors ويا الماتلاب بس في نفس الوقت الغرافز يشتغلون… هذي اول مشكلة
تعبت يا ربي من كثر ما احاول وافشل
والمشكلة العويصة المصيبة الثانية…ان مطلوب مني احول الماتلاب كود الى C++ بس انا في حياتي ما جربت لغة السي بلاس, فواااايد صعب علي اني اسويه ولا حتى عندي فكره عنه… والمصيبة الاكبر يا شباب ان الدكتور مسافر ومب راااجع والله العظيم اني خااايف ساعدوني كيف اقدر احول الماتلاب كودز للسي بلاس؟؟؟؟
وانا بخلي الكودز اللي سويتهم للماتلاب وبليز ساعدوني يا اهل الخير…
feel free to send me private messages or ask me to contact me on email
والله اللي بيساعدني ما بنساله المعروف طول عمري
اول كود مخزن باسم
dif1d.m
[FONT=Courier New][SIZE=2][COLOR=#228b22]%This code solve the one-dimensional heat diffusion equation
%for the problem of a bar which is initially at T-Tinit and
%suddenly the temperatures at the left and right change to
%Tleft and Tright.
%
%Upon discretization in space by a finite difference method,
%the result is a system of ODE's of the form
%
%u_t=Au+b
%
%The code calculates A and b. Then, these are used in
%various different time intergration methods to evolve the
%temperature in time. See lecture notes for more info.
%
%Using this code, we will talk about the following issues:
%
% * Use of explicit methods for stiff systems
% * Implementation of implicit methods for linear systems
% * Sparse vs. full matrices
% * Behavior of trapezoidal when lambda dt -> -infinity.
[/color][/size][/font][FONT=Courier New][SIZE=2]clear [/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]all[/color][/size][/font][FONT=Courier New][SIZE=2];
close [/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]all[/color][/size][/font][FONT=Courier New][SIZE=2];
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Set non-dimensional thermal coefficient
[/color][/size][/font][FONT=Courier New][SIZE=2]k=1.0; [/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% this is really k/(rho*cp)
% Set length of bar
[/color][/size][/font][FONT=Courier New][SIZE=2]L=1.0; [/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]%non-dimensional
% Set initial temperature
[/color][/size][/font][FONT=Courier New][SIZE=2]Tinit=400;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Set left and right temperatures for t>0
[/color][/size][/font][FONT=Courier New][SIZE=2]Tleft=800;
Tright=1000;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Set up grid size
[/color][/size][/font][FONT=Courier New][SIZE=2]Nx=1000;
h=L/Nx;
x=linspace(0,L,Nx+1);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Set timestep size
[/color][/size][/font][FONT=Courier New][SIZE=2]dt=1e-1;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Calculate number of iterations (Nmax) needed to iterate to t=Tmax
[/color][/size][/font][FONT=Courier New][SIZE=2]Tmax=0.5;
Nmax=floor(Tmax/dt);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Initialize a sparse matrix to hold stiffness & identity matrix
[/color][/size][/font][FONT=Courier New][SIZE=2]A=spalloc (Nx-1,Nx-1,3*(Nx-1));
I=speye(Nx-1);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% These commands allocate full (non-sparse) matrices
%A=zeros(Nx-1,Nx-1);
%I=eye(Nx-1);
% Calculate stiffness matrix
[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]for[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#000000] ii=1:Nx-1, [/color]
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]if[/color][/size][/font][FONT=Courier New][SIZE=2] (ii>1),
A(ii,ii-1)=k/h^2;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]if[/color][/size][/font][FONT=Courier New][SIZE=2] (ii<Nx-1),
A(ii,ii+1)=k/h^2;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][SIZE=2]
A(ii,ii)=-2*k/h^2;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Set forcing vector
[/color][/size][/font][FONT=Courier New][SIZE=2]b=zeros (Nx-1,1);
b(1)=k*Tleft/h^2;
b(Nx-1)=k*Tright/h^2;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Set initial vector
[/color][/size][/font][FONT=Courier New][SIZE=2]v=Tinit*ones (Nx-1,1);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]%Iterate in time
[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]for[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#000000] n=1: Nmax,[/color]
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]%Backward Euler
[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]%v=(I-dt*A)\(v+dt*b);
[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]%Forward Euler
[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]%v=v+dt*(A*v+b);
[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]%Trapezoidal
[/color][/size][/font][FONT=Courier New][SIZE=2]rhs=v+0.5*dt*A*v+dt*b;
G=I-0.5*dt*A;
v=G\rhs;
T=[Tleft; v; Tright];
plot (x,T,[/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'*'[/color][/size][/font][FONT=Courier New][SIZE=2]);
axis([0,1,400,1200]);
drawnow;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font]
The second code is saved as
dif1d_fun.m
[FONT=Courier New][SIZE=2][COLOR=#228b22]
% This routine returns the forcing term for
% a one-dimensional heat diffusion problem
% that has been discretized by finite differences.
% Note that the matrix A and the vector b are pre-computed
% in the main driver routine, dif1d_main.m, and passed
% to this function. Then, this function simply returns
% f(v) = A*v+b. So, in reality, this function is
% not specific to 1-d diffusion.
[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]function[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#000000] [f]=dif1d_fun(t, v, A, b)[/color]
f=A*v+b;
[/size][/font]
The third code is saved as
dif1d_main.m
[FONT=Courier New][SIZE=2][COLOR=#228b22]% Matlab script: dif1d_main.m
%
% This code solve the one-dimensional heat diffusion equation
% for the problem of a bar which is initially at T=Tinit and
% suddenly the temperatures at the left and right change to
% Tleft and Tright.
%
% Upon discretization in space by a finite difference method,
% the result is a system of ODE's of the form,
%
% u_t=Au+b
%
% The code calculate A and b. Then, uses one of Matlab's
% ODE inegrators, either ode45 (which is based on a Runge-Kutta
% method and is not designed for stiff problems) or ode15s (which
% is based on an implicit method and is designed for stiff problems).
%
[/color][/size][/font][FONT=Courier New][SIZE=2]clear [/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]all[/color][/size][/font][FONT=Courier New][SIZE=2];
close [/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]all[/color][/size][/font][FONT=Courier New][SIZE=2];
sflag = input([/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'Use stiff integrator? (1=yes, [default=no]):'[/color][/size][/font][FONT=Courier New][SIZE=2]);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Set non-dimensional thermal coefficient
[/color][/size][/font][FONT=Courier New][SIZE=2]k=1.0; [/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% this is really k/(rho*cp)
% Set length of bar
[/color][/size][/font][FONT=Courier New][SIZE=2]L=1.0; [/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% non-dimensional
% Set initial temperature
[/color][/size][/font][FONT=Courier New][SIZE=2]Tinit=400;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Set left right temperatures for t>0
[/color][/size][/font][FONT=Courier New][SIZE=2]Tleft=800;
Tright=1000;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Set up grid size
[/color][/size][/font][FONT=Courier New][SIZE=2]Nx = input([[/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'Enter number of divisions in x-direction: [default='[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]...[/color][/size][/font][FONT=Courier New][SIZE=2][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'51] '[/color][/size][/font][FONT=Courier New][SIZE=2] ]);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]if[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#000000] (isempty(Nx)),[/color]
Nx=51;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][SIZE=2]h=L/Nx;
x=linspace (0,L,Nx+1);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Calculate number of iterations (Nmax) needed to iterate to t=Tmax
[/color][/size][/font][FONT=Courier New][SIZE=2]Tmax=0.5;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Initialize a sparse matrix to hold stiffness & identity matrix
[/color][/size][/font][FONT=Courier New][SIZE=2]A= spalloc(Nx-1,Nx-1,3*(Nx-1));
I= speye(Nx-1);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Calculate stiffness matrix
[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]for[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#000000] ii=1:Nx-1,[/color]
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]if[/color][/size][/font][FONT=Courier New][SIZE=2](ii>1),
A(ii,ii-1)=k/h^2;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]if[/color][/size][/font][FONT=Courier New][SIZE=2](ii<Nx-1),
A(ii,ii+1)=k/h^2;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][SIZE=2]
A(ii,ii)=-2*k/h^2;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Set forcing vector
[/color][/size][/font][FONT=Courier New][SIZE=2]b=zeros(Nx-1,1);
b(1)=k*Tleft/h^2;
b(Nx-1)=k*Tright/h^2;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Set initial vector
[/color][/size][/font][FONT=Courier New][SIZE=2]v0=Tinit*ones(1,Nx-1);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]if[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#000000] (sflag==1),[/color]
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Call ODE15s
[/color][/size][/font][FONT=Courier New][SIZE=2]options = odeset([/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'Jacobian'[/color][/size][/font][FONT=Courier New][SIZE=2],A);
[t,v]=ode15s (@dif1d_fun, [0 Tmax], v0, options, A, b);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]else
[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Call ODE45
[/color][/size][/font][FONT=Courier New][SIZE=2][t,v]=ode45 (@dif1d_fun, [0 Tmax], v0, [], A,b);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Get midpoint value of T and plot vs. time
[/color][/size][/font][FONT=Courier New][SIZE=2]Tmid = v(:,floor (Nx/2));
plot(t,Tmid);
xlabel([/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'t'[/color][/size][/font][FONT=Courier New][SIZE=2]);
ylabel([/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'T at midpoint'[/color][/size][/font][FONT=Courier New][SIZE=2]);
[/size][/font]
The fourth code is saved as
eig_dif1d.m
[FONT=Courier New][SIZE=2][COLOR=#228b22]% Set non-dimensional thermal coefficient
[/color][/size][/font][FONT=Courier New][SIZE=2]k=1.0; [/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% this is really k/(rho*cp)
% Set length of bar
[/color][/size][/font][FONT=Courier New][SIZE=2]L=1.0; [/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% non-dimensional
% Set up grid size
[/color][/size][/font][FONT=Courier New][SIZE=2]Nx=1000;
h=L/Nx;
x=linspace(0,L,Nx+1);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Initialize a sparse matrix to hold stiffness matrix
[/color][/size][/font][FONT=Courier New][SIZE=2]A=spalloc(Nx-1,Nx-1,3*(Nx-1));
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Calculate stiffnesz matrix
[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]for[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#000000] ii=1:Nx-1,[/color]
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]if[/color][/size][/font][FONT=Courier New][SIZE=2](ii>1),
A(ii,ii-1)=k/h^2;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]if[/color][/size][/font][FONT=Courier New][SIZE=2](ii<Nx-1),
A(ii,ii+1)=k/h^2;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][SIZE=2]
A(ii,ii)=-2*k/h^2;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Find eigenvalues of A
[/color][/size][/font][FONT=Courier New][SIZE=2]lam=eig(A);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Determine minimum and maximum eigenvalues of A
[/color][/size][/font][FONT=Courier New][SIZE=2][lammin, imin]= min(abs(lam)); [/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Find minimum magnitude eigenvalue
[/color][/size][/font][FONT=Courier New][SIZE=2][lammax, imax]= max(abs(lam)); [/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Find maximum magnitude eigenvalue
% Determine spectral condition number
[/color][/size][/font][FONT=Courier New][SIZE=2]scn=lammax/lammin;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Print out info to screen
[/color][/size][/font][FONT=Courier New][SIZE=2]fprintf([/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'Minimum eigenvalue = %f
'[/color][/size][/font][FONT=Courier New][SIZE=2],lam(imin));
fprintf([/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'Maximum eigenvalue = %f
'[/color][/size][/font][FONT=Courier New][SIZE=2],lam(imax));
fprintf([/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'Spectral condition number = %f
'[/color][/size][/font][FONT=Courier New][SIZE=2],scn);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Plot eigenvalues
[/color][/size][/font][FONT=Courier New][SIZE=2]plot (real(lam), imag(lam),[/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'*'[/color][/size][/font][FONT=Courier New][SIZE=2]);
xlabel ([/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'Real\lambda'[/color][/size][/font][FONT=Courier New][SIZE=2]);
ylabel ([/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'Imag\lambda'[/color][/size][/font][FONT=Courier New][SIZE=2]);
[/size][/font]
The fifth code is saved as
mstepstab.m
[FONT=Courier New][SIZE=2][COLOR=#228b22]% Specify theta
[/color][/size][/font][FONT=Courier New][SIZE=2]theta=linspace(0,2*pi,301);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Set growth factor
[/color][/size][/font][FONT=Courier New][SIZE=2]g=exp(i*theta);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Calculate lambda dt location for each method and plot.
% 1st order
[/color][/size][/font][FONT=Courier New][SIZE=2]z1=(g-1)./g;
plot(z1);hold [/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]on[/color][/size][/font][FONT=Courier New][SIZE=2];
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Second order
[/color][/size][/font][FONT=Courier New][SIZE=2]z2=(g.^2-4/3*g+1/3)./(2/3*g.^2);
plot(real(z2), imag (z2),[/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'r'[/color][/size][/font][FONT=Courier New][SIZE=2]);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Third order
[/color][/size][/font][FONT=Courier New][SIZE=2]z3=(g.^3 - 18/11*g.^2 +9/11*g - 2/11)./(6/11*g.^3);
plot(z3);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Fourth order
[/color][/size][/font][FONT=Courier New][SIZE=2]z4=(g.^4-48/25*g.^3+36/25*g.^2-16/25*g+3/25)./(12/25*g.^4);
plot(z4);
axis([/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'equal'[/color][/size][/font][FONT=Courier New][SIZE=2]);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]%axis ([-3,3,-3,3]);
[/color][/size][/font][FONT=Courier New][SIZE=2]xlabel ([/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'Real\lambda\Delta t'[/color][/size][/font][FONT=Courier New][SIZE=2]);
ylabel ([/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'Imag\lambda\Delta t'[/color][/size][/font][FONT=Courier New][SIZE=2]);
grid [/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]on[/color][/size][/font][FONT=Courier New][SIZE=2];
[/size][/font]
The sixth file is saved as
stiff_err.m
[FONT=Courier New][SIZE=2]Nvec=[100,500,1000,2000,4000,8000,12000,100000];
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]for[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#000000] nnn=1:length(Nvec),[/color]
clear [/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]ufe[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]ube[/color][/size][/font][FONT=Courier New][SIZE=2];
N=Nvec(nnn)
stiff_forced;
ube= ue - vb;
ebe(nnn)=sqrt (Tmax/N*ube*ube');
ufe= ue - vf;
efe(nnn)=sqrt(Tmax/N*ufe*ufe');
ute= ue - vt;
ete(nnn)=sqrt(Tmax/N*ute*ute');
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][SIZE=2]hvec = Tmax./Nvec;
loglog(hvec,ebe,[/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'*-'[/color][/size][/font][FONT=Courier New][SIZE=2]); hold [/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]on[/color][/size][/font][FONT=Courier New][SIZE=2];
loglog(hvec,efe,[/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'*--'[/color][/size][/font][FONT=Courier New][SIZE=2]);
loglog(hvec,ete,[/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]'*-.'[/color][/size][/font][FONT=Courier New][SIZE=2]);
grid [/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]on[/color][/size][/font][FONT=Courier New][SIZE=2];
[/size][/font]
The seventh and the last code is saved as
stiff_forced.m
[FONT=Courier New][SIZE=2]
clear [/size][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]vb[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]vf[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]vt[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#a020f0]t[/color][/size][/font][FONT=Courier New][SIZE=2];
lambda = -1000;
Tmax = 10;
dt = Tmax/N;
t=linspace(0,Tmax,N+1);
vf=zeros(1,N+1); [/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% It is important to allocate the space for
[/color][/size][/font][FONT=Courier New][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% speed as N gets large.
[/color][/size][/font][FONT=Courier New][SIZE=2]
vb=zeros(1,N+1);
vt=zeros(1,N+1);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Exact solution
[/color][/size][/font][FONT=Courier New][SIZE=2]omega=1;
c2=100/(omega^2+lambda^2);
c1=1+omega*c2;
ue=c1*exp(lambda*t)+c2*(-lambda*sin(omega*t)-omega*cos(omega*t));
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Forward Euler method
[/color][/size][/font][FONT=Courier New][SIZE=2]vf(1)=1.0;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]for[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#000000] n = 1:N,[/color]
vf(n+1)=vf(n)+dt*(lambda*vf(n)+100*sin(t(n)));
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Backward Euler method
[/color][/size][/font][FONT=Courier New][SIZE=2]vb(1)=1.0;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]for[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#000000] n=1:N,[/color]
vb(n+1)=(vb(n)+dt*100*sin(t(n+1)))/(1-dt*lambda);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% Trapezoidal method
[/color][/size][/font][FONT=Courier New][SIZE=2]vt(1)=1.0;
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]for[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#000000] n=1:N,[/color]
vt(n+1)=(vt(n)+0.5*dt*lambda*vt(n)+0.5*dt*100*(sin(t(n))+sin(t(n+1))))/(1-0.5*dt*lambda);
[/size][/font][FONT=Courier New][SIZE=2][COLOR=#0000ff]end
[/color][/size][/font][FONT=Courier New][SIZE=2][COLOR=#228b22]% plot (t,ue,'r');
% hold on;
% plot (t,vt);
% axis ([0,10,-1,1]);
% hold off;
[/color][/size][/font]
people plz note tht the codes written above are all saved in the same directory because they’re related to each other
اسمحولي يا اخوان لان ادري الموضوع شوي طويل بس انا عجزت…لو كان بيدي شي جام ما قصرت
تحياتي لكم وانتظر ردودكم
في امان الله